{ "cells": [ { "cell_type": "code", "execution_count": null, "id": "4761e952-3831-4409-b75f-97a0c1414041", "metadata": { "id": "4761e952-3831-4409-b75f-97a0c1414041" }, "outputs": [], "source": [ "import copy\n", "from types import SimpleNamespace\n", "\n", "import vrplib\n", "import matplotlib\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import numpy.random as rnd\n", "\n", "from alns import ALNS\n", "from alns.accept import RecordToRecordTravel\n", "from alns.select import RouletteWheel\n", "from alns.stop import MaxIterations" ] }, { "cell_type": "code", "execution_count": null, "id": "c8958300-e8ac-42bc-9543-d7ebcb0e364b", "metadata": { "id": "c8958300-e8ac-42bc-9543-d7ebcb0e364b" }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": null, "id": "343632a7", "metadata": { "id": "343632a7" }, "outputs": [], "source": [ "SEED = 1234" ] }, { "cell_type": "markdown", "id": "420d8c99-2ad1-4ec8-a29b-061bce69be91", "metadata": { "id": "420d8c99-2ad1-4ec8-a29b-061bce69be91" }, "source": [ "# The capacitated vehicle routing problem\n", "The vehicle routing problem (VRP) is one of the most studied problems in operations research. Given a fleet of vehicles and customers, the goal is to construct a set of routes such that each customer is visited exactly once, while minimizing the total distance traveled by the vehicles. \n", "\n", "Despite decades of research, the VRP (and variants thereof) remains a very hard problem to solve and new algorithms continue to be developed to address this problem. A related and interesting fact is that ALNS was originally proposed by [Ropke and Pisinger (2006)](https://pubsonline.informs.org/doi/abs/10.1287/trsc.1050.0135?casa_token=-DeLGU-Nr_4AAAAA:hTPxhhAn8TRi5h8p5LdQ_r-tQ1j4lCD4-K4ZR4gSi0e9O6reL6vcyfC0NZmkW1hoQGkUEjcumwH6) to solve many variants of the vehicle routing problem. \n", "\n", "In this notebook, we use ALNS to solve the most famous VRP variant: the *Capacitated Vehicle Routing Problem (CVRP)*. The CVRP can be described using an undirected graph $G=(V,E)$, where $V$ is the vertex set and $E$ is the edge set. The vertex set $V$ is partitioned into $V=\\{0\\} \\cup V_c$, where $0$ is the depot and $V_c=\\{1, \\dots, n\\}$ denotes the set of $n$ customers. Each customer $i \\in V_c$ has a demand $q_i > 0$. A distance $d_{ij}$ is associated with each edge $(i, j) \\in E$. We assume that we have an unlimited fleet of homogeneous vehicles with capacity $Q$ located at the depot. A feasible solution to the CVRP is a set of routes, each served by a single vehicle, such that each customer is served exactly once and none of the routes exceed the vehicle capacity. The goal is to minimize the total distance traveled." ] }, { "cell_type": "markdown", "id": "6ff79ded", "metadata": { "id": "6ff79ded" }, "source": [ "## Data\n", "[CVRPLIB](http://vrp.atd-lab.inf.puc-rio.br/index.php/en/) contains a large collection of CVRP benchmark instances. The library is actively maintained and new best known solutions are updated regularly. We use the `vrplib` package to read the `ORTEC-n242-k12` instance, which consists of 241 customers (+ 1 depot) and 12 vehicles, but we assume that an unlimited number of vehicles is available." ] }, { "cell_type": "code", "execution_count": null, "id": "da0dfed7", "metadata": { "id": "da0dfed7" }, "outputs": [], "source": [ "data = vrplib.read_instance(\"data/ORTEC-n242-k12.vrp\")\n", "bks = SimpleNamespace(**vrplib.read_solution(\"data/ORTEC-n242-k12.sol\"))" ] }, { "cell_type": "markdown", "id": "2398da3e", "metadata": { "id": "2398da3e" }, "source": [ "The `bks` variable contains the best known solution. Let's plot it, together with the coordinates of the customers:" ] }, { "cell_type": "code", "execution_count": null, "id": "776916f5", "metadata": { "id": "776916f5" }, "outputs": [], "source": [ "def plot_solution(solution, name=\"CVRP solution\"):\n", " \"\"\"\n", " Plot the routes of the passed-in solution.\n", " \"\"\"\n", " fig, ax = plt.subplots(figsize=(12, 10))\n", " cmap = matplotlib.cm.rainbow(np.linspace(0, 1, len(solution.routes)))\n", "\n", " for idx, route in enumerate(solution.routes):\n", " ax.plot(\n", " [data[\"node_coord\"][loc][0] for loc in [0] + route + [0]],\n", " [data[\"node_coord\"][loc][1] for loc in [0] + route + [0]],\n", " color=cmap[idx],\n", " marker=\".\",\n", " )\n", "\n", " # Plot the depot\n", " kwargs = dict(label=\"Depot\", zorder=3, marker=\"*\", s=750)\n", " ax.scatter(*data[\"node_coord\"][0], c=\"tab:red\", **kwargs)\n", "\n", " ax.set_title(f\"{name}\\n Total distance: {solution.cost}\")\n", " ax.set_xlabel(\"X-coordinate\")\n", " ax.set_ylabel(\"Y-coordinate\")\n", " ax.legend(frameon=False, ncol=3)" ] }, { "cell_type": "code", "execution_count": null, "id": "f951a07d-ff46-44a3-ba74-844858bfe843", "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 637 }, "id": "8f19d325", "outputId": "2c9bb974-4e8e-4cd1-fad5-668aa93c18fd" }, "outputs": [], "source": [ "plot_solution(bks, name=\"Best known solution\")" ] }, { "cell_type": "markdown", "id": "49a823ca", "metadata": { "id": "49a823ca" }, "source": [ "## Solution state" ] }, { "cell_type": "code", "execution_count": null, "id": "ed998203-e6ee-4c5b-9e52-fe7d99ac3c59", "metadata": { "id": "ed998203-e6ee-4c5b-9e52-fe7d99ac3c59" }, "outputs": [], "source": [ "class CvrpState:\n", " \"\"\"\n", " Solution state for CVRP. It has two data members, routes and unassigned.\n", " Routes is a list of list of integers, where each inner list corresponds to\n", " a single route denoting the sequence of customers to be visited. A route\n", " does not contain the start and end depot. Unassigned is a list of integers,\n", " each integer representing an unassigned customer.\n", " \"\"\"\n", "\n", " def __init__(self, routes, unassigned=None):\n", " self.routes = routes\n", " self.unassigned = unassigned if unassigned is not None else []\n", "\n", " def copy(self):\n", " return CvrpState(copy.deepcopy(self.routes), self.unassigned.copy())\n", "\n", " def objective(self):\n", " \"\"\"\n", " Computes the total route costs.\n", " \"\"\"\n", " return sum(route_cost(route) for route in self.routes)\n", "\n", " @property\n", " def cost(self):\n", " \"\"\"\n", " Alias for objective method. Used for plotting.\n", " \"\"\"\n", " return self.objective()\n", "\n", " def find_route(self, customer):\n", " \"\"\"\n", " Return the route that contains the passed-in customer.\n", " \"\"\"\n", " for route in self.routes:\n", " if customer in route:\n", " return route\n", "\n", " raise ValueError(f\"Solution does not contain customer {customer}.\")\n", "\n", "\n", "def route_cost(route):\n", " distances = data[\"edge_weight\"]\n", " tour = [0] + route + [0]\n", "\n", " return sum(\n", " distances[tour[idx]][tour[idx + 1]] for idx in range(len(tour) - 1)\n", " )" ] }, { "cell_type": "markdown", "id": "06375ce5", "metadata": { "id": "06375ce5" }, "source": [ "## Destroy operators" ] }, { "cell_type": "markdown", "id": "21b6e774-d6f5-4350-bbae-a1358f3e7d18", "metadata": { "id": "21b6e774-d6f5-4350-bbae-a1358f3e7d18" }, "source": [ "Destroy operators break parts of a solution down, leaving an incomplete state. This is the first part of each iteration of the ALNS meta-heuristic; the incomplete solution is subsequently repaired by any one repair operator. We will consider one destroy operator: **random removal**. We will also use a separate parameter, the degree of destruction, to control the extent of the damage done to a solution in each step." ] }, { "cell_type": "code", "execution_count": null, "id": "e34acdfb-c269-46cd-875d-348c7fd45716", "metadata": { "id": "e34acdfb-c269-46cd-875d-348c7fd45716" }, "outputs": [], "source": [ "degree_of_destruction = 0.05\n", "customers_to_remove = int((data[\"dimension\"] - 1) * degree_of_destruction)" ] }, { "cell_type": "code", "execution_count": null, "id": "7c3d8c4a-2a4b-4866-a281-454f7c3bb61c", "metadata": { "id": "7c3d8c4a-2a4b-4866-a281-454f7c3bb61c" }, "outputs": [], "source": [ "def random_removal(state, rng):\n", " \"\"\"\n", " Removes a number of randomly selected customers from the passed-in solution.\n", " \"\"\"\n", " destroyed = state.copy()\n", "\n", " for customer in rng.choice(\n", " range(1, data[\"dimension\"]), customers_to_remove, replace=False\n", " ):\n", " destroyed.unassigned.append(customer)\n", " route = destroyed.find_route(customer)\n", " route.remove(customer)\n", "\n", " return remove_empty_routes(destroyed)\n", "\n", "\n", "def remove_empty_routes(state):\n", " \"\"\"\n", " Remove empty routes after applying the destroy operator.\n", " \"\"\"\n", " state.routes = [route for route in state.routes if len(route) != 0]\n", " return state" ] }, { "cell_type": "markdown", "id": "eb0a7653", "metadata": { "id": "eb0a7653" }, "source": [ "## Repair operators\n", "We implement a simple, **greedy repair** strategy. It iterates over the set of unassigned customers and finds the best route and index to insert to, i.e., with the least increase in cost." ] }, { "cell_type": "code", "execution_count": null, "id": "95537882", "metadata": { "id": "95537882" }, "outputs": [], "source": [ "def greedy_repair(state, rng):\n", " \"\"\"\n", " Inserts the unassigned customers in the best route. If there are no\n", " feasible insertions, then a new route is created.\n", " \"\"\"\n", " rng.shuffle(state.unassigned)\n", "\n", " while len(state.unassigned) != 0:\n", " customer = state.unassigned.pop()\n", " route, idx = best_insert(customer, state)\n", "\n", " if route is not None:\n", " route.insert(idx, customer)\n", " else:\n", " state.routes.append([customer])\n", "\n", " return state\n", "\n", "\n", "def best_insert(customer, state):\n", " \"\"\"\n", " Finds the best feasible route and insertion idx for the customer.\n", " Return (None, None) if no feasible route insertions are found.\n", " \"\"\"\n", " best_cost, best_route, best_idx = None, None, None\n", "\n", " for route in state.routes:\n", " for idx in range(len(route) + 1):\n", "\n", " if can_insert(customer, route):\n", " cost = insert_cost(customer, route, idx)\n", "\n", " if best_cost is None or cost < best_cost:\n", " best_cost, best_route, best_idx = cost, route, idx\n", "\n", " return best_route, best_idx\n", "\n", "\n", "def can_insert(customer, route):\n", " \"\"\"\n", " Checks if inserting customer does not exceed vehicle capacity.\n", " \"\"\"\n", " total = data[\"demand\"][route].sum() + data[\"demand\"][customer]\n", " return total <= data[\"capacity\"]\n", "\n", "\n", "def insert_cost(customer, route, idx):\n", " \"\"\"\n", " Computes the insertion cost for inserting customer in route at idx.\n", " \"\"\"\n", " dist = data[\"edge_weight\"]\n", " pred = 0 if idx == 0 else route[idx - 1]\n", " succ = 0 if idx == len(route) else route[idx]\n", "\n", " # Increase in cost of adding customer, minus cost of removing old edge\n", " return dist[pred][customer] + dist[customer][succ] - dist[pred][succ]" ] }, { "cell_type": "markdown", "id": "63b03225", "metadata": { "id": "63b03225" }, "source": [ "## Initial solution\n", "We need an initial solution that is going to be destroyed and repaired by the ALNS heuristic. To this end, we use a simple *nearest neighbor (NN)* heuristic. NN starts with an empty solution and iteratively adds the nearest customer to the routes. If there are no routes available, then a new route is created." ] }, { "cell_type": "code", "execution_count": null, "id": "56191127", "metadata": { "id": "56191127" }, "outputs": [], "source": [ "def neighbors(customer):\n", " \"\"\"\n", " Return the nearest neighbors of the customer, excluding the depot.\n", " \"\"\"\n", " locations = np.argsort(data[\"edge_weight\"][customer])\n", " return locations[locations != 0]\n", "\n", "\n", "def nearest_neighbor():\n", " \"\"\"\n", " Build a solution by iteratively constructing routes, where the nearest\n", " customer is added until the route has met the vehicle capacity limit.\n", " \"\"\"\n", " routes = []\n", " unvisited = set(range(1, data[\"dimension\"]))\n", "\n", " while unvisited:\n", " route = [0] # Start at the depot\n", " route_demands = 0\n", "\n", " while unvisited:\n", " # Add the nearest unvisited customer to the route till max capacity\n", " current = route[-1]\n", " nearest = [nb for nb in neighbors(current) if nb in unvisited][0]\n", "\n", " if route_demands + data[\"demand\"][nearest] > data[\"capacity\"]:\n", " break\n", "\n", " route.append(nearest)\n", " unvisited.remove(nearest)\n", " route_demands += data[\"demand\"][nearest]\n", "\n", " customers = route[1:] # Remove the depot\n", " routes.append(customers)\n", "\n", " return CvrpState(routes)" ] }, { "cell_type": "code", "execution_count": null, "id": "1ceb1d6e", "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 637 }, "id": "1ceb1d6e", "outputId": "173c9006-baa1-4c96-bfae-9deb731212cd" }, "outputs": [], "source": [ "plot_solution(nearest_neighbor(), \"Nearest neighbor solution\")" ] }, { "cell_type": "markdown", "id": "ff1bc37f", "metadata": { "id": "ff1bc37f" }, "source": [ "## Heuristic solution\n", "\n", "Let's now construct our ALNS heuristic. Since we only have one destroy and repair operator, we do not actually use any adaptive operator selection -- but you can easily add more destroy and repair operators. " ] }, { "cell_type": "code", "execution_count": null, "id": "f850b236", "metadata": { "id": "f850b236" }, "outputs": [], "source": [ "alns = ALNS(rnd.default_rng(SEED))\n", "\n", "alns.add_destroy_operator(random_removal)\n", "\n", "alns.add_repair_operator(greedy_repair)" ] }, { "cell_type": "code", "execution_count": null, "id": "c9da2b44", "metadata": { "id": "c9da2b44" }, "outputs": [], "source": [ "num_iterations = 3000\n", "init = nearest_neighbor()\n", "select = RouletteWheel([25, 5, 1, 0], 0.8, 1, 1)\n", "accept = RecordToRecordTravel.autofit(\n", " init.objective(), 0.02, 0, num_iterations\n", ")\n", "stop = MaxIterations(num_iterations)\n", "\n", "result = alns.iterate(init, select, accept, stop)" ] }, { "cell_type": "code", "execution_count": null, "id": "d06ff996", "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 674 }, "id": "d06ff996", "outputId": "019c9c13-f487-4c67-adc7-597c3fbdbe7c" }, "outputs": [], "source": [ "solution = result.best_state\n", "objective = solution.objective()\n", "pct_diff = 100 * (objective - bks.cost) / bks.cost\n", "\n", "print(f\"Best heuristic objective is {objective}.\")\n", "print(\n", " f\"This is {pct_diff:.1f}% worse than the optimal solution, which is {bks.cost}.\"\n", ")" ] }, { "cell_type": "code", "execution_count": null, "id": "f309ec78-2dfa-4d5e-9d9d-7355357ed382", "metadata": {}, "outputs": [], "source": [ "_, ax = plt.subplots(figsize=(12, 6))\n", "result.plot_objectives(ax=ax)" ] }, { "cell_type": "markdown", "id": "cfe7b5a8-983f-4dce-b538-fd846c0520f9", "metadata": {}, "source": [ "Let's have a look at the solution:" ] }, { "cell_type": "code", "execution_count": null, "id": "97a55721-4627-48bc-bddd-3eee492e28f7", "metadata": {}, "outputs": [], "source": [ "plot_solution(solution, \"Simple ALNS\")" ] }, { "cell_type": "markdown", "id": "P8O54LTGjjBu", "metadata": { "id": "P8O54LTGjjBu" }, "source": [ "## Slack-induced substring removal\n", "The simple destroy and repair operator from above work fine, but there are better destroy and repair operators for CVRP. One example is the *Slack Induction by Substring Removal (SISR)* method proposed by [Christiaens and Vanden Berghe (2020)](https://pubsonline.informs.org/doi/abs/10.1287/trsc.2019.0914?casa_token=lPUXU1Ax8PIAAAAA:yTE9Pu6L9QGPRu_vt-ZMHF0AZvL9gV0fNS4QAUTOJboQcgTVyOR9_RTbm9rZcImyKI4GUW9pLv1j). SISR obtains state-of-the-art results using a destroy operator that, instead of removing random customers, removes partial routes (called *strings*) that are all located near each other. Moreover, a blinking feature is added to the greedy repair operator, where certain insertion checks are skipped. For more details, we refer to the paper.\n", "\n", "In the following, we will implement a simplified version of the string removal operator and replace the random destroy operator." ] }, { "cell_type": "code", "execution_count": null, "id": "59de74b2-7dbb-4893-ba1a-57f895685d69", "metadata": { "id": "WlBmF77PIDGy" }, "outputs": [], "source": [ "MAX_STRING_REMOVALS = 2\n", "MAX_STRING_SIZE = 12\n", "\n", "\n", "def string_removal(state, rng):\n", " \"\"\"\n", " Remove partial routes around a randomly chosen customer.\n", " \"\"\"\n", " destroyed = state.copy()\n", "\n", " avg_route_size = int(np.mean([len(route) for route in state.routes]))\n", " max_string_size = max(MAX_STRING_SIZE, avg_route_size)\n", " max_string_removals = min(len(state.routes), MAX_STRING_REMOVALS)\n", "\n", " destroyed_routes = []\n", " center = rng.integers(1, data[\"dimension\"])\n", "\n", " for customer in neighbors(center):\n", " if len(destroyed_routes) >= max_string_removals:\n", " break\n", "\n", " if customer in destroyed.unassigned:\n", " continue\n", "\n", " route = destroyed.find_route(customer)\n", " if route in destroyed_routes:\n", " continue\n", "\n", " customers = remove_string(route, customer, max_string_size, rng)\n", " destroyed.unassigned.extend(customers)\n", " destroyed_routes.append(route)\n", "\n", " return destroyed\n", "\n", "\n", "def remove_string(route, cust, max_string_size, rng):\n", " \"\"\"\n", " Remove a string that constains the passed-in customer.\n", " \"\"\"\n", " # Find consecutive indices to remove that contain the customer\n", " size = rng.integers(1, min(len(route), max_string_size) + 1)\n", " start = route.index(cust) - rng.integers(size)\n", " idcs = [idx % len(route) for idx in range(start, start + size)]\n", "\n", " # Remove indices in descending order\n", " removed_customers = []\n", " for idx in sorted(idcs, reverse=True):\n", " removed_customers.append(route.pop(idx))\n", "\n", " return removed_customers" ] }, { "cell_type": "code", "execution_count": null, "id": "SeJWYINg20Fn", "metadata": { "id": "SeJWYINg20Fn" }, "outputs": [], "source": [ "alns = ALNS(rnd.default_rng(SEED))\n", "\n", "alns.add_destroy_operator(string_removal)\n", "\n", "alns.add_repair_operator(greedy_repair)" ] }, { "cell_type": "code", "execution_count": null, "id": "Rx2-Q0V-28Lx", "metadata": { "id": "Rx2-Q0V-28Lx" }, "outputs": [], "source": [ "num_iterations = 3000\n", "init = nearest_neighbor()\n", "select = RouletteWheel([25, 5, 1, 0], 0.8, 1, 1)\n", "accept = RecordToRecordTravel.autofit(\n", " init.objective(), 0.02, 0, num_iterations\n", ")\n", "stop = MaxIterations(num_iterations)\n", "\n", "result = alns.iterate(init, select, accept, stop)" ] }, { "cell_type": "code", "execution_count": null, "id": "V8oOY04m3mf3", "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 404 }, "id": "V8oOY04m3mf3", "outputId": "9d222fce-d1c1-44f2-dbc3-b3f6a4636532" }, "outputs": [], "source": [ "_, ax = plt.subplots(figsize=(12, 6))\n", "result.plot_objectives(ax=ax)" ] }, { "cell_type": "code", "execution_count": null, "id": "XLrbh5loClBZ", "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 674 }, "id": "XLrbh5loClBZ", "outputId": "1e3b5ca8-49d1-4ad6-fcc3-449958839953" }, "outputs": [], "source": [ "solution = result.best_state\n", "objective = solution.objective()\n", "pct_diff = 100 * (objective - bks.cost) / bks.cost\n", "\n", "print(f\"Best heuristic objective is {objective}.\")\n", "print(\n", " f\"This is {pct_diff:.1f}% worse than the optimal solution, which is {bks.cost}.\"\n", ")\n", "plot_solution(solution, \"String removals\")" ] }, { "cell_type": "markdown", "id": "d541a8e7-461a-40ab-9d85-2c94712c1c29", "metadata": {}, "source": [ "## Conclusions" ] }, { "cell_type": "markdown", "id": "7424ca3f-b4e7-459d-a307-be2d26426cdc", "metadata": {}, "source": [ "In this notebook we implemented two heuristics for the CVRP, using the ALNS meta-heuristic framework. The first heuristic used a random customer destroy operator and we obtained a solution which is 10% worse than the best known solution. The second heuristic used a destroy operator which removes strings arround a randomly selected customer and we obtained a solution that is only 3% worse than the best known solution. \n", "\n", "This example shows that by constructing problem-specific operators, one can create even more powerful ALNS heuristics." ] } ], "metadata": { "colab": { "collapsed_sections": [ "68965045-48f2-41ac-ae4e-394ae450c5de" ], "provenance": [] }, "kernelspec": { "display_name": ".venv", "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.9.19" } }, "nbformat": 4, "nbformat_minor": 5 }