Using AequilibraE and Google OR-Tools to compute the TSP

TSP with AequilibraE.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from os.path import join\n",
    "from aequilibrae.project import Project\n",
    "from aequilibrae.matrix import AequilibraeData, AequilibraeMatrix\n",
    "from aequilibrae.paths import NetworkSkimming, SkimResults, PathResults\n",
    "import numpy as np\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ortools.constraint_solver import routing_enums_pb2\n",
    "from ortools.constraint_solver import pywrapcp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Computing a cost matrix between stops"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The example project comes from one of the standard AequilibraE example projects \n",
    "# https://www.aequilibrae.com/data/Chicago.7z\n",
    "\n",
    "fldr = 'D:/net_tests/'\n",
    "proj_name = 'chicago_imported.sqlite'\n",
    "\n",
    "project = Project()\n",
    "project.load(join(fldr, proj_name))\n",
    "project.network.build_graphs()\n",
    "\n",
    " # we grab the graph for cars\n",
    "graph = project.network.graphs['c']\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Getting 30 nodes at RANDOM here, but this is up to you!\n",
    "targets = 30\n",
    "\n",
    "curr = project.conn.cursor()\n",
    "curr.execute('select node_id from nodes')\n",
    "all_nodes = [a[0] for a in curr.fetchall()]\n",
    "\n",
    "nodes_for_skims = np.random.choice(all_nodes, targets)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We prepare the graph for skimming\n",
    "graph.prepare_graph(nodes_for_skims)\n",
    "graph.set_graph('free_flow_time') # let's say we want to minimize time\n",
    "graph.set_skimming(['free_flow_time']) # And will skim time and distance\n",
    "\n",
    "# We are using regular nodes, so the links that feed into the centroids should be preserved\n",
    "graph.set_blocked_centroid_flows(False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run skimming\n",
    "res = SkimResults()\n",
    "res.prepare(graph)\n",
    "\n",
    "res.compute_skims()\n",
    "skm = res.skims\n",
    "mat = skm.get_matrix('free_flow_time')\n",
    "mat = (mat*100).astype(int) # set the matrix to integers, which is a requirement from ortools"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Travelling-Salesman Problem solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# little help function to retrieve the route using the network IDs\n",
    "def get_route(manager, routing):\n",
    "    index = routing.Start(0)\n",
    "    nodes = []\n",
    "    while not routing.IsEnd(index):\n",
    "        p = skm.index[manager.IndexToNode(index)]\n",
    "        nodes.append(p)\n",
    "        previous_index = index\n",
    "        index = solution.Value(routing.NextVar(index))\n",
    "    \n",
    "    nodes.append(skm.index[manager.IndexToNode(index)])\n",
    "    return nodes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The actual TSP computation\n",
    "\n",
    "depot = 0 # It is the matrix index that counts\n",
    "vehicles = 1\n",
    "\n",
    "# Create the routing index manager.\n",
    "manager = pywrapcp.RoutingIndexManager(mat.shape[0], vehicles, depot)\n",
    "\n",
    "# Create Routing Model.\n",
    "routing = pywrapcp.RoutingModel(manager)\n",
    "\n",
    "def distance_callback(from_index, to_index):\n",
    "    # Convert from routing variable Index to distance matrix NodeIndex.\n",
    "    from_node = manager.IndexToNode(from_index)\n",
    "    to_node = manager.IndexToNode(to_index)\n",
    "    \n",
    "    return mat[from_node, to_node]\n",
    "\n",
    "transit_callback_index = routing.RegisterTransitCallback(distance_callback)\n",
    "\n",
    "# Define cost of each arc.\n",
    "routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)\n",
    "\n",
    "# Setting first solution heuristic.\n",
    "search_parameters = pywrapcp.DefaultRoutingSearchParameters()\n",
    "search_parameters.first_solution_strategy = (routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC)\n",
    "\n",
    "# Solve the problem.\n",
    "solution = routing.SolveWithParameters(search_parameters)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Retrieve the solution\n",
    "if solution:\n",
    "    n = get_route(manager, routing)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Re-building the path through the network"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_links = []\n",
    "pthres = PathResults()\n",
    "pthres.prepare(graph)\n",
    "for i in range(len(n)-1):\n",
    "    f = n[i]\n",
    "    t = n[i + 1]\n",
    "    pthres.compute_path(f, t)\n",
    "    all_links.extend(list(pthres.path))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import geopandas as gpd\n",
    "import folium"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sql = \"SELECT link_id, a_node as in_route, Hex(ST_AsBinary(GEOMETRY)) as geom FROM links;\"\n",
    "route = gpd.GeoDataFrame.from_postgis(sql, project.conn,  geom_col=\"geom\")\n",
    "\n",
    "# We make a field to identify those in route\n",
    "route.loc[:, 'in_route'] = 0\n",
    "route.loc[route.link_id.isin(all_links),'in_route'] = 1.0\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "order = pd.DataFrame({'node_id':n[:-1],\n",
    "                       'sequence':np.arange(1,len(n))})\n",
    "\n",
    "sql = \"SELECT node_id, Hex(ST_AsBinary(GEOMETRY)) as geom FROM nodes;\"\n",
    "node_layer = gpd.GeoDataFrame.from_postgis(sql, project.conn,  geom_col=\"geom\")\n",
    "stops = node_layer[node_layer.node_id.isin(n)]\n",
    "stops = stops.merge(order, on='node_id')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "routelines = route.to_json()\n",
    "routepoints = stops.to_json()\n",
    "\n",
    "node_ids = folium.features.GeoJsonTooltip(['node_id','sequence'])\n",
    "\n",
    "\n",
    "colors = {0:'blue', 1:'red'}\n",
    "def line_style_function(feature):\n",
    "    return {'opacity': 1, \n",
    "            'weight': 3 * float(feature['properties']['in_route']) + 0.1, \n",
    "            'color': colors[int(feature['properties']['in_route'])]}\n",
    "\n",
    "\n",
    "mapa = folium.Map([41.883718, -87.632382],\n",
    "                  zoom_start=8,\n",
    "                  tiles='cartodbpositron')\n",
    "\n",
    "lines = folium.features.GeoJson(routelines, name='Route', style_function=line_style_function).add_to(mapa)\n",
    "points = folium.features.GeoJson(routepoints, name='Stops', tooltip=node_ids).add_to(mapa)\n",
    "\n",
    "mapa"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}

This is what the resulting map looks like