{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Exercises in Transportation Optimization\n", "---\n", "Niels Lindner, Ricardo Euler, Pedro Maristany de las Casas, Stephan Schwartz
\n", "CO@Work 2020 - September 22, 2020\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Line Planning\n", "---\n", "\n", "### Input:\n", "\n", "* an undirected graph $G = (V, E)$, vertices are stations of a public transportation network,\n", "* an OD matrix $(d_{st}) \\in \\mathbb R^{V \\times V}$ estimating the number of passengers traveling from $s$ to $t$,\n", "* a line pool $\\mathcal L_0$ of paths (lines) in $G$,\n", "* a capacity $C \\geq 0$, maximum number of passengers per vehicle,\n", "* a frequency bound $F \\in \\mathbb N_0$ on the maximum number of vehicles per edge,\n", "* travel time $t \\in \\mathbb R_{\\geq 0}^E$, passengers' travel time on each edge,\n", "* operational line costs $c \\in \\mathbb R_{\\geq 0}^E$, costs for operating a vehicle on an edge,\n", "* a scalarization parameter $\\lambda \\in \\mathbb R_{\\geq 0}$.\n", "\n", "### Output:\n", "A line plan $(\\mathcal L, f)$, consisting of a subset $\\mathcal L \\subseteq \\mathcal L_0$ of lines, and a frequency $f_\\ell \\in \\mathbb N_0$ for each $\\ell \\in \\mathcal L$ such that\n", "* the vehicles of the lines in $\\mathcal L$ cover the passenger demand,\n", "* the frequencies $f$ satisfy the frequency bound on each edge,\n", "* passenger travel time + $\\lambda \\cdot$ operational line costs is minimized." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Task 0: Resolving data inconsistencies\n", "---\n", "As in real world data, there are some errors in our data. Please replace lines/lines-star.csv and lines/lines-2020.csv with the updated versions from here:\n", "https://cloud.zib.de/index.php/s/G29QG6ZAojGDQMg\n", "You can do this by going back to the directory view, and clicking \"Upload\" on the top right corner." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Task 1: Reading the input graph\n", "---\n", "Create an instance of a networkx.Graph from the stations in the file lines/stations.csv and the edges in lines/lines-links.csv. The resulting undirected graph should have 23 vertices and 30 edges. Store the x and y fiels for later use.\n", " " ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "vertices: 23, edges: 30\n" ] } ], "source": [ "import csv\n", "import networkx\n", "\n", "graph = networkx.Graph()\n", "\n", "with open(\"lines/stations.csv\") as f:\n", " reader = csv.DictReader(f, delimiter=\";\")\n", " for row in reader:\n", " graph.add_node(row[\"id\"], x=float(row[\"x\"]), y=float(row[\"y\"]))\n", " \n", "with open(\"lines/lines-links.csv\") as f:\n", " reader = csv.DictReader(f, delimiter=\";\")\n", " for row in reader:\n", " graph.add_edge(*row[\"stations\"].split(\",\"))\n", " \n", "print(f\"vertices: {len(graph.nodes)}, edges: {len(graph.edges)}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Task 2: Reading the OD matrix\n", "---\n", "Read the OD information from od.csv and store it in a double dict." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "od = {s: {t: 0 for t in graph.nodes} for s in graph.nodes}\n", "\n", "with open(\"lines/od.csv\") as f:\n", " reader = csv.DictReader(f, delimiter=\";\")\n", " for row in reader:\n", " od[row[\"source\"]][row[\"target\"]] = float(row[\"demand\"])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Task 3: Reading the line pool\n", "---\n", "Construct a line pool from the lines in lines/lines-star.csv. The pool should be stored as a dict whose keys are from the id field in the csv, and each value is a list corresponding to the edge sequence of the line.\n", " " ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "{'1': [('Gvc', 'Ut'), ('Ut', 'Apd'), ('Apd', 'Hgl'), ('Hgl', 'Odz')], '2': [('Gv', 'Shl'), ('Shl', 'Asdz'), ('Asdz', 'Ut'), ('Ut', 'Ah'), ('Ah', 'Zv')], '3': [('Asd', 'Ut'), ('Ut', 'Ehv'), ('Ehv', 'Std'), ('Std', 'Mt')], '4': [('Rsd', 'Rtd'), ('Rtd', 'Bd'), ('Bd', 'Ut'), ('Ut', 'Zl'), ('Zl', 'Hr'), ('Hr', 'Lw')], '5': [('Lls', 'Zl'), ('Zl', 'Asn'), ('Asn', 'Gn')]}\n" ] } ], "source": [ "pool = {}\n", "\n", "with open(\"lines/lines-star.csv\") as f:\n", " reader = csv.DictReader(f, delimiter=\";\")\n", " for row in reader:\n", " path = row[\"stations\"].split(\",\")\n", " pool[row[\"id\"]] = list(zip(path, path[1:]))\n", " \n", "print(pool)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Task 4: Setting up the MIP model\n", "---\n", "\n", "We model the line planning problem as follows:\n", "\n", "We turn $G$ into a directed graph $G' = (V, A)$, replacing each undirected edge $\\{v, w\\}$ by two antiparallel arcs $vw$ and $wv$. For an arc $vw \\in A$, set $\\mathcal L_{vw} := \\mathcal L_{\\{v,w\\}} := \\{\\ell \\in \\mathcal L_0: \\{v, w\\} \\in \\ell\\}$. We further extend $t$ and $c$ to $A$, keeping the value of the undirected edges. Finally define $D := \\{(s,t) \\mid d_{st} > 0\\}$.\n", "\n", "\\begin{align}\n", "&\\text{Minimize} & \\sum_{a \\in A} \\sum_{(s, t) \\in D} t_a x_{st, a} + \\lambda \\sum_{a \\in A} \\sum_{\\ell \\in \\mathcal L_a} c_a f_\\ell&\\\\\n", "&\\text{s.t.} & \\sum_{(s, t) \\in D} x_{st,a} - C \\cdot \\sum_{\\ell \\in \\mathcal L_a} f_\\ell &\\leq 0, &a \\in A,\\\\\n", "& & \\sum_{a \\in \\delta^+(v)} x_{st,a} - \\sum_{a \\in \\delta^-(v)} x_{st,a} &= 0, &(s, t) \\in D, v \\in V \\setminus \\{s, t\\},\\\\\n", "& & \\sum_{a \\in \\delta^+(s)} x_{st,a} - \\sum_{a \\in \\delta^-(s)} x_{st,a} &= d_{st}, &(s, t) \\in D,\\\\\n", "& & \\sum_{a \\in \\delta^+(t)} x_{st,a} - \\sum_{a \\in \\delta^-(t)} x_{st,a} &= -d_{st}, &(s, t) \\in D,\\\\\n", "& & \\sum_{\\ell \\in \\mathcal L_e} f_\\ell &\\leq F, &e \\in E,\\\\\n", "& & x_{st,a} &\\geq 0, &(s, t) \\in \\mathcal D, a \\in A,\\\\\n", "& & f_\\ell &\\geq 0, &\\ell \\in \\mathcal L_0,\\\\\n", "& & f_\\ell &\\in \\mathbb Z, &\\ell \\in \\mathcal L_0.\n", "\\end{align}\n", "\n", "The routing of the passengers is hence done on the directed version $G'$ of $G$, while lines stay undirected.\n", "\n", "Task: Set up this MIP model for the above input graph, OD matrix, and line pool. We set $C := 200$, $F := 10$, $\\lambda := 100$. For an edge $e = \\{v,w\\} \\in E$, we set $t_e := c_e := \\sqrt{(x_v-x_w)^2+(y_v-y_w)^2}$. You may use pyscipopt or alternatively write directly a .lp-File in CPLEX LP file format." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:4: UserWarning: linked SCIP 7.0 is not recommended for this version of PySCIPOpt - use version 7.0.1\n", " after removing the cwd from sys.path.\n" ] } ], "source": [ "#initialize pyscipopt\n", "import pyscipopt\n", "\n", "model = pyscipopt.Model()" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "#input parameters and helper functions\n", "C = 200\n", "F = 10\n", "lam = 100\n", "\n", "import itertools\n", "od_pairs = [(s, t) for (s, t) in itertools.permutations(graph.nodes, 2) if od[s][t] > 0]\n", "\n", "import math\n", "def cost(graph, v, w):\n", " return math.sqrt(pow(graph.nodes[v][\"x\"] - graph.nodes[w][\"x\"], 2) \\\n", " + pow(graph.nodes[v][\"y\"] - graph.nodes[w][\"y\"], 2))\n", "\n", "def lines_on_edge(graph, v, w):\n", " return [l for l in pool if (v, w) in pool[l] or (w, v) in pool[l]]" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "#variables\n", "x = {}\n", "f = {}\n", "\n", "for (s, t) in od_pairs:\n", " for (v, w) in graph.edges:\n", " x[s,t,v,w] = model.addVar(f\"x_{s},{t},{v},{w}\", vtype=\"C\", lb=0, ub=None, obj=cost(graph, v, w))\n", " x[s,t,w,v] = model.addVar(f\"x_{s},{t},{w},{v}\", vtype=\"C\", lb=0, ub=None, obj=cost(graph, v, w))\n", "\n", "for l in pool:\n", " f[l] = model.addVar(f\"f_{l}\", vtype=\"I\", lb=0, ub=None, \\\n", " obj=lam*sum([cost(graph, v, w) for (v, w) in pool[l]]))" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "wrote problem to file b'model.lp'\n" ] } ], "source": [ "#flow-frequency coupling constraints\n", "for (v, w) in graph.edges:\n", " model.addCons(pyscipopt.quicksum(x[s,t,v,w] for (s, t) in od_pairs) \\\n", " <= C * pyscipopt.quicksum(f[l] for l in lines_on_edge(graph, v, w)))\n", " model.addCons(pyscipopt.quicksum(x[s,t,w,v] for (s, t) in od_pairs) \\\n", " <= C * pyscipopt.quicksum(f[l] for l in lines_on_edge(graph, v, w)))\n", "\n", "#flow conservation constraints\n", "for (s, t) in od_pairs:\n", " for v in graph.nodes:\n", " b = 0\n", " if v == s:\n", " b = od[s][t]\n", " elif v == t:\n", " b = -od[s][t]\n", " model.addCons(pyscipopt.quicksum(x[s,t,v,w] for w in graph.neighbors(v)) \\\n", " - pyscipopt.quicksum(x[s,t,w,v] for w in graph.neighbors(v)) == b)\n", "\n", "#frequency bound constraints \n", "for (v, w) in graph.edges:\n", " model.addCons(pyscipopt.quicksum(f[l] for l in lines_on_edge(graph, v, w)) <= F) \n", " \n", "#write to file\n", "model.writeProblem(\"model.lp\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Task 5: Solving the MIP\n", "---\n", "\n", "Invoke SCIP and let it solve the problem for you. What is the total cost of the line plan? Print the line frequencies and the total number of passengers on each undirected edge." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "presolving:\n", "(round 1, fast) 4582 del vars, 2100 del conss, 0 add conss, 2757 chg bounds, 0 chg sides, 0 chg coeffs, 0 upgd conss, 0 impls, 0 clqs\n", "(round 2, fast) 6791 del vars, 3373 del conss, 0 add conss, 7067 chg bounds, 12 chg sides, 0 chg coeffs, 0 upgd conss, 0 impls, 0 clqs\n", "(round 3, fast) 8982 del vars, 3703 del conss, 0 add conss, 7623 chg bounds, 26 chg sides, 0 chg coeffs, 0 upgd conss, 0 impls, 0 clqs\n", "(round 4, fast) 9474 del vars, 3720 del conss, 0 add conss, 7629 chg bounds, 40 chg sides, 0 chg coeffs, 0 upgd conss, 0 impls, 0 clqs\n", "(round 5, fast) 9483 del vars, 3724 del conss, 0 add conss, 7629 chg bounds, 44 chg sides, 0 chg coeffs, 0 upgd conss, 0 impls, 0 clqs\n", "presolving (6 rounds: 6 fast, 1 medium, 1 exhaustive):\n", " 9485 deleted vars, 3724 deleted constraints, 0 added constraints, 7629 tightened bounds, 0 added holes, 44 changed sides, 0 changed coefficients\n", " 0 implications, 0 cliques\n", "presolved problem has 0 variables (0 bin, 0 int, 0 impl, 0 cont) and 0 constraints\n", "Presolving Time: 0.09\n", "\n", " time | node | left |LP iter|LP it/n|mem/heur|mdpt |vars |cons |rows |cuts |sepa|confs|strbr| dualbound | primalbound | gap | compl. \n", "t 0.1s| 1 | 0 | 0 | - | trivial| 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 2.901111e+06 | 2.901111e+06 | 0.00%| unknown\n", "\n", "SCIP Status : problem is solved [optimal solution found]\n", "Solving Time (sec) : 0.11\n", "Solving Nodes : 1\n", "Primal Bound : +2.90111137641328e+06 (1 solutions)\n", "Dual Bound : +2.90111137641328e+06\n", "Gap : 0.00 %\n" ] } ], "source": [ "#solve\n", "model.redirectOutput()\n", "model.optimize()" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Total cost of line plan: 2901111.3764132783\n" ] } ], "source": [ "sol = model.getBestSol()\n", "\n", "#total cost of line plan\n", "print(f\"Total cost of line plan: {model.getObjVal()}\")" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1: 4.0\n", "2: 2.0\n", "3: 4.0\n", "4: 7.0\n", "5: 3.0\n" ] } ], "source": [ "#line frequencies\n", "for l in pool:\n", " print(f\"{l}: {sol[f[l]]}\")" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Ah -> Ut : 310.0\n", "Ah -> Zv : 206.0\n", "Apd -> Asd : 0.0\n", "Apd -> Hgl : 266.0\n", "Apd -> Ut : 627.0\n", "Asd -> Lls : 0.0\n", "Asd -> Shl : 0.0\n", "Asd -> Ut : 37.0\n", "Asdz -> Lls : 0.0\n", "Asdz -> Shl : 163.0\n", "Asdz -> Ut : 306.0\n", "Asn -> Gn : 203.0\n", "Asn -> Zl : 485.0\n", "Bd -> Rtd : 269.0\n", "Bd -> Ut : 648.0\n", "Ehv -> Std : 412.0\n", "Ehv -> Ut : 652.0\n", "Gv -> Gvc : 0.0\n", "Gv -> Rtd : 0.0\n", "Gv -> Shl : 123.0\n", "Gvc -> Shl : 0.0\n", "Gvc -> Ut : 241.0\n", "Hgl -> Odz : 47.0\n", "Hr -> Lw : 153.0\n", "Hr -> Zl : 336.0\n", "Lls -> Zl : 123.0\n", "Mt -> Std : 239.0\n", "Rsd -> Rtd : 112.0\n", "Rtd -> Ut : 0.0\n", "Ut -> Zl : 1258.0\n" ] } ], "source": [ "#total flow per edge\n", "for (v, w) in graph.edges:\n", " print(f\"{v:4s} -> {w:4s}: {sum([sol[x[s,t,v,w]] for (s, t) in od_pairs]):6.1f}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The routing is rather boring on a tree. There are, e.g., 1258 passengers traveling fom Ut to Zl. The MIP should become infeasible if the total capacity on this edge is below 1258, e.g., if infrastructure still allows $F = 10$ trains per hour, but shorter vehicles with capacity $C = 125$ are desired. For $C = 126$, the problem becomes feasible again, but the cost of the line plan is higher (why?). If you want, you can test this out and play around a little." ] }, { "cell_type": "code", "execution_count": 91, "metadata": {}, "outputs": [], "source": [ "model.freeTransform()\n", "model.freeProb()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Task 6: A larger line pool\n", "---\n", "We consider the same network and parameters, but another line pool, given in lines/lines-2020.csv. What is different in the solution process? Is the obtained line plan cheaper?\n", " " ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "#with open(\"lines/lines-2020.csv\") as f:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Finale\n", "---\n", "An advanced, but state-of-the-art method is to use column generation for determining the lines, so that one is not restricted to a chosen line pool. The pricing problem is in general a constrained shortest path problem. In an optimal solution, only few lines will have positive frequencies, and this is why column generation can be fruitful. For bigger instances, writing down the full model can even be infeasible.\n", "\n", "Our instance is small enough to avoid column generation by generating all simple paths for the line pool. This can be done with networkx.all_simple_paths(graph, source, target)." ] }, { "cell_type": "code", "execution_count": 100, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "7302\n" ] } ], "source": [ "import itertools\n", "\n", "pool = {}\n", "i = 0\n", "\n", "for (v, w) in itertools.combinations(graph.nodes, 2): \n", " for path in networkx.all_simple_paths(graph, v, w):\n", " i += 1\n", " pool[i] = list(zip(path, path[1:]))\n", " \n", "print(len(pool))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Observe that the resulting MIP model is already really large, making it hard to solve for SCIP (solving takes more than one hour). In fact, this is overkill, because we can reach the same result with $\\mathcal L_0 = E$, reducing the computation time to a few seconds." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "30\n" ] } ], "source": [ "import itertools\n", "\n", "pool = {}\n", "i = 0\n", "\n", "for (v, w) in graph.edges:\n", " i += 1\n", " pool[i] = [(v, w)]\n", " \n", "print(len(pool))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " The line plan has again become cheaper. If you write the resulting lines into a CSV, you can draw the lines with the lines/draw.py script." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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.6.9" } }, "nbformat": 4, "nbformat_minor": 4 }