Solving traveling salesman problems using Solver Foundation

Update: see the comments below for some helpful hints. If you are unable to run this with your version of Solver Foundation and Gurobi, consider installing the lp_solve plugin for MSF. More details on this thread.

Here’s an example that I walked through during yesterday’s INFORMS session.  Erwin has two blog postings about Solver Foundation and the traveling salesman problem, but I want to throw in my two cents because I want to emphasize a couple of points:

  1. By combining C# and Solver Foundation Services it is possible to express complex models clearly and succinctly.
  2. It is very easy to build powerful, reusable model libraries using C# and Solver Foundation Services.
  3. Solver Foundation Services code can be used in many different application environments (ASP.Net, silverlight, DB, command line apps, WPF, …) with minimal changes.

The traveling salesman problem is a classical problem in computer science, and you should bow your head in shame if you don’t know about it (and turn in your conference badge if you happen to be in Phoenix). William Cook’s book on the traveling salesman problem is a wonderful read. A salesperson needs to make a tour of a number of cities.  The restrictions are that she wants to visit each city once and only once, and she wants to minimize the distance travelled.  This is perhaps the definitive example of an NP-hard problem.

TSP can be solved using mixed integer programming – optimizing a linear goal with linear constraints, where some of the decision variables are integer.  In this first post I will show how to formulate and solve a TSP model using Solver Foundation Services.  In my second post I will show how to use the Gurobi MIP solver using SFS.   There are many different ways to model the TSP – here is a nice introduction.  My goal is to provide a clear, complete example – not build a “production level” TSP model, so I am going to choose a model formulation that dates back to 1960!  First, I need to establish a couple of building blocks that will help me construct the data for the model.  We need to know the distances between each pair of cities.  Typically we are provided the coordinates of the cities and need to derive the distances.  So I will introduce a Coordinate class that contains properties for the (x, y) coordinates, and properties to convert to latitude and longitude.  Finally, a method that computes the distance between points.

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using Microsoft.SolverFoundation.Services;

namespace Microsoft.SolverFoundation.Samples {
  class TravelingSalesman {
    // TSP coordinate.
    public class Coordinate {
      public int Name { get; set; }

      // X-coordinate (from TSPLIB)
      public double X { get; set; }

      // Y-coordinate (from TSPLIB)
      public double Y { get; set; }

      public Coordinate(int name, double x, double y) {
        Name = name;
        X = x;
        Y = y;
      }

      // Latitude in radians.
      public double Latitude {
        get { return Math.PI * (Math.Truncate(X) + 5 * (X - Math.Truncate(X)) / 3) / 180; }
      }

      // Longitude in radians.
      public double Longitude {
        get { return Math.PI * (Math.Truncate(Y) + 5 * (Y - Math.Truncate(Y)) / 3) / 180; }
      }

      // Geographic distance between two points (as an integer).
      public int Distance(Coordinate p) {
        double q1 = Math.Cos(Longitude - p.Longitude);
        double q2 = Math.Cos(Latitude - p.Latitude);
        double q3 = Math.Cos(Latitude + p.Latitude);
        // There may rounding difficulties her if the points are close together...just sayin'.
        return (int)(6378.388 * Math.Acos(0.5 * ((1 + q1) * q2 - (1 - q1) * q3)) + 1);
      }
    }

    // TSP city-city arc.
    public class Arc {
      public int City1 { get; set; }
      public int City2 { get; set; }
      public double Distance { get; set; }
    }

    // Burma14 from TSPLIB. Optimal tour = 3323.
    private static Coordinate[] data = new Coordinate[] {
      new Coordinate(0, 16.47, 96.10),
      new Coordinate(1, 16.47, 94.44),
      new Coordinate(2, 20.09, 92.54),
      new Coordinate(3, 22.39, 93.37),
      new Coordinate(4, 25.23, 97.24),
      new Coordinate(5, 22.00, 96.05),
      new Coordinate(6, 20.47, 97.02),
      new Coordinate(7, 17.20, 96.29),
      new Coordinate(8, 16.30, 97.38),
      new Coordinate(9, 14.05, 98.12),
      new Coordinate(10, 16.53, 97.38),
      new Coordinate(11, 21.52, 95.59),
      new Coordinate(12, 19.41, 97.13),
      new Coordinate(13, 20.09, 94.55)
    };
}

(The data for this 14-city problem comes from the TSPLIB library). If you’ve been following my blog you know that the building blocks of a Solver Foundation model are: sets, parameters, decisions, goals, and constraints. I am going to implement a simple formulation that is centered around the following (indexed) decisions:

  • Assign[i,j]: this is equal to 1 if the optimal tour contains a trip (or arc) from city i to city j.
  • Rank[i]: this is equal to the number of cities visited after arriving at city i.

We have one parameter in our model:

  • Distance[I,j]: the distance from city i to city j.

With that in mind, here’s the model.  Explanation of the goals and constraints follow.

    
public static void Run() {
      SolverContext context = SolverContext.GetContext();
      context.ClearModel();
      Model model = context.CreateModel();

      // ------------
      // Parameters
      Set city = new Set(Domain.IntegerNonnegative, "city");
      Parameter dist = new Parameter(Domain.Real, "dist", city, city);
      var arcs = from p1 in data
                 from p2 in data
                 select new Arc { City1 = p1.Name, City2 = p2.Name, Distance = p1.Distance(p2) };
      dist.SetBinding(arcs, "Distance", "City1", "City2");
      model.AddParameters(dist);

      // ------------
      // Decisions
      Decision assign = new Decision(Domain.IntegerRange(0, 1), "assign", city, city);
      Decision rank = new Decision(Domain.RealNonnegative, "rank", city);
      model.AddDecisions(assign, rank);

      // ------------
      // Goal: minimize the length of the tour.
      Goal goal = model.AddGoal("TourLength", GoalKind.Minimize,
        Model.Sum(Model.ForEach(city, i => Model.ForEachWhere(city, j => dist[i, j] * assign[i, j], j => i != j))));

      // ------------
      // Enter and leave each city only once.
      int N = data.Length;
      model.AddConstraint("assign1",
        Model.ForEach(city, i => Model.Sum(Model.ForEachWhere(city, j => assign[i, j],
          j => i != j)) == 1));
      model.AddConstraint("assign2",
        Model.ForEach(city, j => Model.Sum(Model.ForEachWhere(city, i => assign[i, j], i => i != j)) == 1));

      model.AddConstraint("A1", Model.ForEach(city, i => Model.Sum(Model.ForEachWhere(city, j => assign[i, j], j => i != j)) == 1));
      model.AddConstraint("A2", Model.ForEach(city, j => Model.Sum(Model.ForEachWhere(city, i => assign[i, j], i => i != j)) == 1));

      // Forbid subtours (Miller, Tucker, Zemlin - 1960...)
      model.AddConstraint("nosubtours",
        Model.ForEach(city,
          i => Model.ForEachWhere(city,
            j => rank[i] + 1 <= rank[j] + N * (1 - assign[i, j]),
            j => Model.And(i != j, i >= 1, j >= 1)
          )
        )
      );

      Solution solution = context.Solve();

      // Retrieve solution information.
      Console.WriteLine("Cost = {0}", goal.ToDouble());
      Console.WriteLine("Tour:");
      var tour = from p in assign.GetValues() where (double)p[0] > 0.9 select p[2];
      foreach (var i in tour.ToArray()) {
        Console.Write(i + " -> ");
      }
      Console.WriteLine();
    }

In my humble opinion, the “Parameter data =” line is an awesome example of the power of LINQ data binding in Solver Foundation.  We generate the 2D matrix of distances using a single LINQ expression. It would be incredibly easy to change the code to retrieve the coordinate data from a database (perhaps using a LINQ expression once again), a file, or even a user application.

The goal is straightforward: minimize the distance traveled.  This is a product of the selected arcs and the distance matrix.   We have two types of constraints:

  • Assignment constraints: these ensure that we enter and leave each city only once.
  • Subtour constraints: these ensure that we do not have any subtours. In a four city problem {A, B, C, D}, for example, we cannot have two cycles (A, B), (C, D). We need to have one tour that contains all the cities.

The assignment constraints are easy using the ForEach and ForEachWhere operations.  I use ForEachWhere because I want to disallow arcs that enter and leave the same city – that doesn’t make sense.  The subtour constraint is a little more complicated. It relates the “assign” and “rank” decisions. The key fact is that if there is an arc from city i to city j, rank[i] + 1 == j. Of course, if the (i, j) arc is not part of the optimal tour then all bets are off. Last note: notice that I can mix parameters, decisions, and C# variables in my expressions.

Getting the cost is very easy using goal.ToDouble().  We can get the tour using either Assign or Rank.  I have chosen to use Assign because it gives me another opportunity to use LINQ.  When you call GetValues() on a decision, you get arrays that contain the value along with the indexes for each decision.  In this case, the last entry in the array is the one we are interested in. There are other ways to conveniently query decsision results, I’ll save that for another time.

The next post will show how we can use Solver Foundation’s plug-in model to tune the behavior of the Gurobi MIP solver.

28 comments

  1. Dear Nathan,
    Thank you very much for the sources and explanation. I’ve got a question, the obtained solution is optimal? Solver Foundation provide exact solution and it is not a heuristic?
    Best regards,
    Sergey

  2. Hi Sergey, thanks for your comment. Yes, the obtained solution is provably optimal. Results returned by the linear and mixed integer linear solvers are globally optimal (unless the solution status says otherwise, for example if the problem is infeasible).

    Nate

    1. Dear Nathan,
      Thank you for the answer.

      The code works perfectly well as it gives optimal result. But the sequence seems not to be correct. E.g. optimal is: 335, but when I summarize the distances between the obtained nodes in sequence, then it is 783.
      I think the optimality is found correct, I tested this code and compared results with Brute-force algorithms, for 6 points. But likely there is an error of extracting the right tour from the solver. It seems like extracted tour is not last found or so.

      I would be very grateful if you could give any suggestions.
      Best regards,
      Sergey

  3. The solution printed to Console does not appear to be correct or is at least misleading. The output:
    5 -> 10 -> 15 -> …
    implies you travel from city5 to city10 to city15…

    But the output code appears to be written as a list of destinations ordered by origin city with no such that:
    city5 comes after 0, city10 comes after 1, city15 comes after 3, etc.

    -Ed

  4. Dear Nathan,
    I used the report command:

    “Report report = solution.GetReport();
    Console.Write(“{0}”, report);”

    And it is possible to extract the right sequence out of there from “assign”.
    Could you recommend any other elegant way?

    Strange, but
    ” var tour = from p in assign.GetValues() where (double)p[0] > 0.9 select p[2];”
    gives something wrong.
    Best regards,
    Sergey

    1. The problem is in the display of the result. This seems to work:

      var tour = from p in assign.GetValues() where (double)p[0] > 0.9 select p;
      foreach (var t in tour)
      {
      Console.WriteLine(“{0} -> {1}”, t[1], t[2]);
      }

      -Ed

  5. I tried to run this for 25 cities, but it exceeded the Gurobi solver throttle. Is there a way to tell SFS to use a different solver? Is there a good FOSS solver out there with a plugin for SFS?

    Thanks,
    -Ed

      1. natebrix :
        You can try the lpsolve connector. This thread has more:
        http://social.msdn.microsoft.com/Forums/en/solverfoundation/thread/b500d84f-05ce-4a87-af51-3f70be51b443

        Thanks you for the answer. I attached Ipsolver into the project.
        I got this from report:
        “Solve Time (ms): 15412
        Total Time (ms): 15593
        Solve Completion Status: Optimal
        Solver Selected: Microsoft.SolverFoundation.Solvers.SimplexSolver”
        This is for 15 cities configuration. Gurobi solved it within 500 ms. IPSolver 30 times slower.

        Ed, how is it in your example? Slow too?

        Best regards,
        Sergey

      2. natebrix :
        You can try the lpsolve connector. This thread has more:
        http://social.msdn.microsoft.com/Forums/en/solverfoundation/thread/b500d84f-05ce-4a87-af51-3f70be51b443

        Somehow I managed to force Gurobi Solver to work with 30 cities.
        May be the last time I did something wrong, or now it works, because I got an academic license.
        The time needed for 30 points is:
        Solve Time (ms): 46261
        Total Time (ms): 46460
        Solve Completion Status: Optimal
        Solver Selected: SolverFoundation.Plugin.Gurobi.GurobiSolver

  6. Hi Sergey,

    I have not yet been able to solve using the LPSolve plugin. If I figure it out, I’ll follow up with performance observations. Right now I’m getting the following:

    “A call to PInvoke function ‘LpSolvePlugIn!LpSolveNativeInterface.lpsolve::set_timeout’ has unbalanced the stack. This is likely because the managed PInvoke signature does not match the unmanaged target signature. Check that the calling convention and parameters of the PInvoke signature match the target unmanaged signature.”

    It is occuring in the plugin code at:
    lpsolve.set_timeout(_lp, prms._LpSolveTimeout);

    Thanks,
    -Ed

    1. I got it going, and I’m seeing a similar performance difference. The Gurobi solved 20 cities in about 2 seconds, and the LP took about 23.

      1. Ed, I would recommend you to get an academic license from Gurobi (in case if it suits you) and try with it. But I did not use DLL from their package, but rather from here: http://blogs.msdn.com/b/solverfoundation/archive/2011/06/17/gurobi-4-5-connector-for-solver-foundation-3-0.aspx
        It seems that this library also start to work, when you got a license.
        In case if it does not work, try these DLLs: https://dl.dropbox.com/u/30020830/gurobi45.dll

        https://dl.dropbox.com/u/30020830/GurobiPlugIn.dll

        And put them here:
        C:\Program Files\Microsoft Solver Foundation\3.0.2.10889\Plugins
        Although I have 64bit PC, this win32 library works fine. Have no idea how does it work.

        Let me know, when you download libs, so that I can delete them from DropBox.

  7. Thanks Sergey,

    I’d prefer to find/use a good open source solver and it turns out that the LPSolve is actually working marginally faster for me than Gurobi in the hacked up test case with which I’m playing. Here are two runs on 20 cities:

    Version: Microsoft Solver Foundation 3.0.2.10889 Express Edition
    Model Name: DefaultModel
    Capabilities Applied: MILP
    Solve Time (ms): 45
    Total Time (ms): 466
    Solve Completion Status: Optimal
    Solver Selected: SolverFoundation.Plugin.LpSolve.LpSolveSolver

    Version: Microsoft Solver Foundation 3.0.2.10889 Express Edition
    Model Name: DefaultModel
    Capabilities Applied: MILP
    Solve Time (ms): 76
    Total Time (ms): 534
    Solve Completion Status: Optimal
    Solver Selected: SolverFoundation.Plugin.Gurobi.GurobiSolver

    I’m going to go back through the sample and try to identify the minimal changes to get the LPSolve performing similar to the Gurboi, but two things that may have helped are changing the distance parameter from Domain.Real to Domain.IntegerNonnegative and adding this directive:
    context.Solve(new SolverFoundation.Plugin.LpSolve.LpSolveDirective());

    -Ed

  8. Hi, Nate, how can one change the code so that it solves the classic TSP, where the salesperson begins and ends in city 0, so that we would have something like 0->3 -> 8 -> … ->0?
    Thanks in advance!

    1. Oh, now I got it. This is the piece of code I was looking for:

      var tour = from p in assign.GetValues() where (double)p[0] > 0.9 select p[2];
      Console.Write(“0 -> “);
      int indice = 0;
      int currentCity;
      for (int i = 0; i < data.Length; i++)
      {
      currentCity = Convert.ToInt32(tour.ToArray()[indice]);
      Console.Write("{0}" + ((i ” :””), currentCity);
      indice = currentCity;
      }
      Console.WriteLine();

      The result is the correct sequence of cities:
      0 -> 1 -> 13 -> 2 -> 3 -> 4 -> 5 -> 11 -> 6 -> 12 -> 7 -> 10 -> 8 -> 9 -> 0

      1. Hi Ricardo,
        But what is the difference, where to start from?
        Imagine the sequences: 0-1-2-3-4-0 and 1-2-3-4-0-1
        I assume that the tour is actualy the same. You can start anywhere you’d like to.
        Am I right?
        Best regards,
        Sergey

      2. Hi, Sergey, you are mathematically correct. But if the salesman lives in city 0, it is only natural for him to start and finish in his own city :)

  9. Ricardo Padua :
    Hi, Sergey, you are mathematically correct. But if the salesman lives in city 0, it is only natural for him to start and finish in his own city

    After you obtain a tour, you can just make some manipulations with the order and assign any starting point you want. :)

  10. Hi Nathan,

    This is excellent code and solved something I truly over complicated. I am using this for solving production line changeover matrixes that are asymmetrical–i.e. if you are making pudding, the line can change over from vanilla to chocolate with a 15 minute wipe down but when going from chocolate to vanilla it takes 8 hours of clean up.

    Anyway, my 20X20 matrix ran great and actually cut 13% off of what I thought was the optimal!

    What do you think the scalability for this is. We are building up a 43X43 matrix now and want to make sure we are capable of doing 100X100.

    Do you or any of the other readers think this will be a problem or have any suggestions before we attempt this?

    Steve

  11. I tried your solution on a Graph with some not existed edges and your code did not provide any useful output. Can u please present a TSP solution for situation like I face?

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s