tsp implementations in or_tools (single core only, ORTools v9.9.3903, Ubuntu 22.04 LTS WSL2)

tsp_ortools.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Comparison of the performance of the Circuit-constraint\n",
    "\n",
    "There are multiple ways to model the Traveling Salesman Problem (TSP) in OR-Tools.\n",
    "In this notebook, we compare the performance of some."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "TIMELIMIT_S = 10\n",
    "STEP_SIZE = 10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate random instances for a benchmark\n",
    "\n",
    "To compare the performance of the different models, we need to generate random instances.\n",
    "Because purely random graphs are not very interesting, we use points in the plane and the distance between them as edge weights.\n",
    "This gives the instances a more realistic structure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import networkx as nx\n",
    "import random\n",
    "import itertools\n",
    "import typing\n",
    "import datetime\n",
    "from ortools.math_opt.python import mathopt, result, callback\n",
    "from ortools.sat.python import cp_model\n",
    "from ortools.sat import sat_parameters_pb2\n",
    "\n",
    "\n",
    "# Generate a random graph with a specified number of nodes.\n",
    "def generate_random_euclidean_graph(\n",
    "    n: int,\n",
    ") -> typing.Tuple[nx.Graph, typing.List[typing.Tuple[int, int]]]:\n",
    "    \"\"\"\n",
    "    Generate a random graph with a specified number of nodes.\n",
    "    The nodes will be integers 0, 1, ..., n-1.\n",
    "    The graph will be a complete graph with Euclidean distance as edge weight.\n",
    "    \"\"\"\n",
    "    points = [(random.randint(0, 10_000), random.randint(0, 10_000)) for _ in range(n)]\n",
    "    G = nx.Graph()\n",
    "    G.add_nodes_from(range(n))\n",
    "    # Complete graph with Euclidean distance as edge weight.\n",
    "    for i, j in itertools.combinations(range(n), 2):\n",
    "        p = points[i]\n",
    "        q = points[j]\n",
    "        dist = round((p[0] - q[0]) ** 2 + (p[1] - q[1]) ** 2)\n",
    "        G.add_edge(i, j, weight=dist)\n",
    "    return G, points"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model 1: AddCircuit-constraint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def solve_with_add_circuit(G: nx.Graph) -> typing.Tuple[bool, float]:\n",
    "    \"\"\"Solve TSP with add_circuit.\"\"\"\n",
    "    model = cp_model.CpModel()\n",
    "\n",
    "    # Variables\n",
    "    edge_vars = dict()\n",
    "    for u, v in G.edges:\n",
    "        edge_vars[u, v] = model.new_bool_var(f\"edge_{u}_{v}\")\n",
    "        edge_vars[v, u] = model.new_bool_var(f\"edge_{v}_{u}\")\n",
    "\n",
    "    # Constraints\n",
    "    # Because the nodes in the graph a indices 0, 1, ..., n-1, we can use the\n",
    "    # indices directly in the constraints. Otherwise, we would have to use\n",
    "    # a mapping from the nodes to indices.\n",
    "    circuit = [(u, v, x) for (u, v), x in edge_vars.items()]\n",
    "    model.add_circuit(circuit)\n",
    "\n",
    "    # Objective\n",
    "    model.minimize(sum(x * G[u][v][\"weight\"] for (u, v), x in edge_vars.items()))\n",
    "\n",
    "    solver = cp_model.CpSolver()\n",
    "    solver.parameters.max_time_in_seconds = TIMELIMIT_S\n",
    "    solver.parameters.linearization_level = 2\n",
    "    solver.parameters.cp_model_presolve = False\n",
    "    solver.parameters.num_search_workers = 1\n",
    "    # solver.parameters.log_search_progress = True\n",
    "    status = solver.solve(model)\n",
    "    return status == cp_model.OPTIMAL, solver.objective_value\n",
    "\n",
    "\n",
    "assert solve_with_add_circuit(generate_random_euclidean_graph(10)[0])[\n",
    "    0\n",
    "], \"Small graph should be solvable.\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model 2: CPSAT Linear Solve\n",
    "(relaxed model, may or may not give optimum results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cpsat_linear_solve(G: nx.Graph) -> typing.Tuple[bool, float]:\n",
    "    model = mathopt.Model(name=\"tsp\")\n",
    "\n",
    "    # Variables\n",
    "    edge_vars = dict()\n",
    "    for u, v in G.edges:\n",
    "        # symmetric variables\n",
    "        x = model.add_binary_variable(name=f\"edge_{u}_{v}\")\n",
    "        edge_vars[u, v] = x\n",
    "        edge_vars[v, u] = x\n",
    "\n",
    "    # Constraints\n",
    "    # Every nodes has two incident edges.\n",
    "    for v in G.nodes:\n",
    "        model.add_linear_constraint(sum(edge_vars[u, v] for u in G.neighbors(v)) == 2)\n",
    "\n",
    "    # Objective\n",
    "    tour_cost = sum(G.edges[u, v][\"weight\"] * edge_vars[(u, v)] for u, v in G.edges)\n",
    "    model.minimize_linear_objective(tour_cost)\n",
    "    cpsat_params = sat_parameters_pb2.SatParameters(num_search_workers=1,cp_model_presolve=False,linearization_level=2)\n",
    "    params = mathopt.SolveParameters(time_limit=datetime.timedelta(seconds=TIMELIMIT_S),cp_sat=cpsat_params)\n",
    "    output = mathopt.solve(model,mathopt.SolverType.CP_SAT,params=params)\n",
    "    objective_value = output.objective_value() if output.has_primal_feasible_solution() else None\n",
    "    return output.termination.reason == result.TerminationReason.OPTIMAL, objective_value\n",
    "\n",
    "\n",
    "assert cpsat_linear_solve(\n",
    "    generate_random_euclidean_graph(20)[0]\n",
    "), \"Small graph should be solvable.\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model 3: SCIP with cuts\n",
    "(relaxed model, may or may not give optimum results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def scip_linear_solve_cuts(G: nx.Graph) -> typing.Tuple[bool, float]:\n",
    "    model = mathopt.Model(name=\"tsp\")\n",
    "\n",
    "    # Variables\n",
    "    edge_vars = dict()\n",
    "    for u, v in G.edges:\n",
    "        # symmetric variables\n",
    "        x = model.add_binary_variable(name=f\"edge_{u}_{v}\")\n",
    "        edge_vars[u, v] = x\n",
    "        edge_vars[v, u] = x\n",
    "\n",
    "    def subtour_callback(cb_data = callback.CallbackData) -> callback.CallbackResult:\n",
    "        cb_res = callback.CallbackResult()\n",
    "        connected_components = list(\n",
    "            nx.connected_components(\n",
    "                nx.Graph(\n",
    "                    (u, v)\n",
    "                    for (u, v), x in edge_vars.items()\n",
    "                    if cb_data.solution[x] == 1.\n",
    "                )\n",
    "            )\n",
    "        )\n",
    "        if len(connected_components) > 1:\n",
    "            for comp in connected_components:\n",
    "                outgoing_edges = sum(\n",
    "                    edge_vars[u, v]\n",
    "                    for u in comp\n",
    "                    for v in G.neighbors(u)\n",
    "                    if v not in comp\n",
    "                )\n",
    "                cb_res.add_generated_constraint(outgoing_edges >= 2,is_lazy=False)\n",
    "        return cb_res\n",
    "\n",
    "    # Constraints\n",
    "    # Every nodes has two incident edges.\n",
    "    for v in G.nodes:\n",
    "        model.add_linear_constraint(sum(edge_vars[u, v] for u in G.neighbors(v)) == 2)\n",
    "\n",
    "    # Objective\n",
    "    tour_cost = sum(G.edges[u, v][\"weight\"] * edge_vars[(u, v)] for u, v in G.edges)\n",
    "    model.minimize_linear_objective(tour_cost)\n",
    "    params = mathopt.SolveParameters(time_limit=datetime.timedelta(seconds=TIMELIMIT_S))\n",
    "    cb_reg = callback.CallbackRegistration()\n",
    "    cb_reg.events.add(callback.Event.MIP_NODE)\n",
    "    cb_reg.add_cuts = True\n",
    "    output = mathopt.solve(model,mathopt.SolverType.GSCIP,params=params,callback_reg=cb_reg,cb=subtour_callback)\n",
    "    objective_value = output.objective_value() if output.has_primal_feasible_solution() else None\n",
    "    return output.termination.reason == result.TerminationReason.OPTIMAL, objective_value\n",
    "\n",
    "\n",
    "assert scip_linear_solve_cuts(\n",
    "    generate_random_euclidean_graph(20)[0]\n",
    "), \"Small graph should be solvable.\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model 4: SCIP with lazy constraints"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def scip_linear_solve_lazy(G: nx.Graph) -> typing.Tuple[bool, float]:\n",
    "    model = mathopt.Model(name=\"tsp\")\n",
    "\n",
    "    # Variables\n",
    "    edge_vars = dict()\n",
    "    for u, v in G.edges:\n",
    "        # symmetric variables\n",
    "        x = model.add_binary_variable(name=f\"edge_{u}_{v}\")\n",
    "        edge_vars[u, v] = x\n",
    "        edge_vars[v, u] = x\n",
    "\n",
    "    def subtour_callback(cb_data = callback.CallbackData) -> callback.CallbackResult:\n",
    "        cb_res = callback.CallbackResult()\n",
    "        connected_components = list(\n",
    "            nx.connected_components(\n",
    "                nx.Graph(\n",
    "                    (u, v)\n",
    "                    for (u, v), x in edge_vars.items()\n",
    "                    if cb_data.solution[x] == 1.\n",
    "                )\n",
    "            )\n",
    "        )\n",
    "        if len(connected_components) > 1:\n",
    "            for comp in connected_components:\n",
    "                outgoing_edges = sum(\n",
    "                    edge_vars[u, v]\n",
    "                    for u in comp\n",
    "                    for v in G.neighbors(u)\n",
    "                    if v not in comp\n",
    "                )\n",
    "                cb_res.add_lazy_constraint(outgoing_edges >= 2)\n",
    "        return cb_res\n",
    "\n",
    "    # Constraints\n",
    "    # Every nodes has two incident edges.\n",
    "    for v in G.nodes:\n",
    "        model.add_linear_constraint(sum(edge_vars[u, v] for u in G.neighbors(v)) == 2)\n",
    "\n",
    "    # Objective\n",
    "    tour_cost = sum(G.edges[u, v][\"weight\"] * edge_vars[(u, v)] for u, v in G.edges)\n",
    "    model.minimize_linear_objective(tour_cost)\n",
    "    params = mathopt.SolveParameters(time_limit=datetime.timedelta(seconds=TIMELIMIT_S))\n",
    "    cb_reg = callback.CallbackRegistration()\n",
    "    cb_reg.events.add(callback.Event.MIP_SOLUTION)\n",
    "    cb_reg.add_lazy_constraints = True\n",
    "    output = mathopt.solve(model,mathopt.SolverType.GSCIP,params=params,callback_reg=cb_reg,cb=subtour_callback)\n",
    "    objective_value = output.objective_value() if output.has_primal_feasible_solution() else None\n",
    "    return output.termination.reason == result.TerminationReason.OPTIMAL, objective_value\n",
    "\n",
    "\n",
    "assert scip_linear_solve_lazy(\n",
    "    generate_random_euclidean_graph(20)[0]\n",
    "), \"Small graph should be solvable.\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model 5: SCIP with no callback\n",
    "(relaxed model, may or may not give optimum results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def scip_linear_solve(G: nx.Graph) -> typing.Tuple[bool, float]:\n",
    "    model = mathopt.Model(name=\"tsp\")\n",
    "\n",
    "    # Variables\n",
    "    edge_vars = dict()\n",
    "    for u, v in G.edges:\n",
    "        # symmetric variables\n",
    "        x = model.add_binary_variable(name=f\"edge_{u}_{v}\")\n",
    "        edge_vars[u, v] = x\n",
    "        edge_vars[v, u] = x\n",
    "\n",
    "    # Constraints\n",
    "    # Every nodes has two incident edges.\n",
    "    for v in G.nodes:\n",
    "        model.add_linear_constraint(sum(edge_vars[u, v] for u in G.neighbors(v)) == 2)\n",
    "\n",
    "    # Objective\n",
    "    tour_cost = sum(G.edges[u, v][\"weight\"] * edge_vars[(u, v)] for u, v in G.edges)\n",
    "    model.minimize_linear_objective(tour_cost)\n",
    "    params = mathopt.SolveParameters(time_limit=datetime.timedelta(seconds=TIMELIMIT_S))\n",
    "    output = mathopt.solve(model,mathopt.SolverType.GSCIP,params=params)\n",
    "    objective_value = output.objective_value() if output.has_primal_feasible_solution() else None\n",
    "    return output.termination.reason == result.TerminationReason.OPTIMAL, objective_value\n",
    "\n",
    "\n",
    "assert scip_linear_solve(\n",
    "    generate_random_euclidean_graph(20)[0]\n",
    "), \"Small graph should be solvable.\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model 6: OR-Tools Routing Library\n",
    "(search strategies which can escape local minima e.g. GUIDED_LOCAL_SEARCH keeps on running even for smaller graphs, if we use GREEDY_DESCENT it stops before reaching optimum even for smaller graphs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ortools.constraint_solver import pywrapcp, routing_enums_pb2\n",
    "\n",
    "def create_distance_callback(G: nx.Graph,manager: pywrapcp.RoutingIndexManager):\n",
    "    def distance_callback(from_index,to_index):\n",
    "        from_node = manager.IndexToNode(from_index)\n",
    "        to_node = manager.IndexToNode(to_index)\n",
    "        return G.edges[from_node,to_node][\"weight\"]\n",
    "    return distance_callback\n",
    "\n",
    "def solve_with_routing_library(G: nx.Graph) -> typing.Tuple[bool, float]:\n",
    "    \"\"\"Solve TSP with add_circuit.\"\"\"\n",
    "    manager = pywrapcp.RoutingIndexManager(G.number_of_nodes(),1,0)\n",
    "    routing = pywrapcp.RoutingModel(manager)\n",
    "    distance_callback = create_distance_callback(G,manager)\n",
    "    transit_callback = routing.RegisterTransitCallback(distance_callback)\n",
    "    routing.SetArcCostEvaluatorOfAllVehicles(transit_callback)\n",
    "    search_parameters = pywrapcp.DefaultRoutingSearchParameters()\n",
    "    search_parameters.local_search_metaheuristic = (routing_enums_pb2.LocalSearchMetaheuristic.GREEDY_DESCENT)\n",
    "    search_parameters.time_limit.FromSeconds(TIMELIMIT_S)\n",
    "    solution = routing.SolveWithParameters(search_parameters)\n",
    "    return routing.solver().WallTime()/1000 < TIMELIMIT_S, solution.ObjectiveValue()\n",
    "\n",
    "\n",
    "assert solve_with_routing_library(generate_random_euclidean_graph(10)[0])[\n",
    "    0\n",
    "], \"Small graph should be solvable.\"\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n=10\n",
      "  add_circuit solved it with objective 73150502.0.\n",
      "  scip_linear_solve_lazy solved it with objective 73150502.0.\n",
      "n=20\n",
      "  add_circuit solved it with objective 98531220.0.\n",
      "  scip_linear_solve_lazy solved it with objective 98531220.00000001.\n",
      "n=30\n",
      "  add_circuit solved it with objective 64131848.0.\n",
      "  scip_linear_solve_lazy solved it with objective 64131848.00000001.\n",
      "n=40\n",
      "  add_circuit solved it with objective 78161598.0.\n",
      "  scip_linear_solve_lazy solved it with objective 78161598.0.\n",
      "n=50\n",
      "  add_circuit solved it with objective 92732234.0.\n",
      "  scip_linear_solve_lazy solved it with objective 92732234.0.\n",
      "n=60\n",
      "  add_circuit solved it with objective 84177646.0.\n",
      "  scip_linear_solve_lazy dropped out.\n",
      "n=70\n",
      "  add_circuit solved it with objective 80133784.0.\n",
      "n=80\n",
      "  add_circuit dropped out.\n"
     ]
    }
   ],
   "source": [
    "# Only allowed accurate models\n",
    "# Disabled routing library as it doesn't give optimal solution\n",
    "challengers = {\n",
    "    \"add_circuit\": solve_with_add_circuit,\n",
    "    # \"cpsat_linear_solve\": cpsat_linear_solve,\n",
    "    # \"scip_linear_solve_cuts\": scip_linear_solve_cuts,\n",
    "    \"scip_linear_solve_lazy\": scip_linear_solve_lazy,\n",
    "    # \"scip_linear_solve\": scip_linear_solve,\n",
    "    # \"solve_with_routing_library\": solve_with_routing_library,\n",
    "}\n",
    "# multiplying initial step by 5 because checking small instances are not necessary\n",
    "n = STEP_SIZE\n",
    "while challengers:\n",
    "    print(f\"n={n}\")\n",
    "    G, points = generate_random_euclidean_graph(n)\n",
    "    dropouts = []\n",
    "    for name, solver in challengers.items():\n",
    "        solved, obj = solver(G)\n",
    "        if not solved:\n",
    "            dropouts.append(name)\n",
    "        else:\n",
    "            print(f\"  {name} solved it with objective {obj}.\")\n",
    "    for name in dropouts:\n",
    "        print(f\"  {name} dropped out.\")\n",
    "        del challengers[name]\n",
    "    n += STEP_SIZE"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "mo310",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}

Hi Arnab, I tried running this. But I got an error. Looks like somehow math_opt is not installed. I'm using python 3.11.0 with ortools==9.8.3296, Windows 10 Enterprise. How can I install ORTools v9.9.3903?

@murali731 You have to build it from source (main branch of ortools git repo). I have a 100MB+ wheel which I can share with you. But I don't want to share such a large file here in the comment. Also, Github doesn't allow python wheels to be uploaded in the comments.