{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "#
Block 2: Linear programming: min-cost flow problems
\n", "###
Alfred Galichon (NYU & Sciences Po)
\n", "##
'math+econ+code' masterclass on optimal transport and economic applications
\n", "
© 2018-2021 by Alfred Galichon. Past and present support from NSF grant DMS-1716489, ERC grant CoG-866274, and contributions by Jules Baudet, Pauline Corblet, Gregory Dannay, and James Nesbit are acknowledged.
\n", "\n", "####
With python code examples
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Learning Objectives \n", "\n", "* Basic concepts of directed networks\n", "\n", "* The min-cost flow problem\n", "\n", "* Duality, optimality and equilibrium\n", "\n", "### References\n", "\n", "* Galichon, *Optimal Transport Methods in Economics*. Ch.8.\n", "\n", "* Tolstoi (1930). Methods of finding the minimal total kilometrage in cargo transportation planning in space. *Transportation Planning* [in Russian]\n", "\n", "* Koopmans (1949). Optimum utilization of the transportation system. *Econometrica*.\n", "\n", "* Schrijver (2002). On the History of the Transportation and Maximum Flow Problems. *Mathematical programming*." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Motivation\n", "\n", "### Motivation: network flow problems\n", "\n", "In 1930 Tolstoi, a Russian engineer, has to optimize the shipping of cement from factories to cities in the Soviet union through railway.\n", "\n", "* each factory produces a fixed number of tons\n", "\n", "* each city needs a fixed number of tons -- for now, we'll assume total production=total consumption\n", "\n", "* each factory is connected by rail with a few cities, and the corresponding distance is given\n", "\n", "* how to ship in order to minimize the total distance travelled?\n", "\n", "This problem belongs to the class of *min-cost flow problems*, an important class of linear programming problems, which are the focus of today's lecture. A decade before the invention of linear programming and the work of Kantorovich, Koopmans and Dantzig, Tolstoi described a heuristic method for solving the problem, which led to the optimal solution. As one can recover using modern tools, and will see that his solution was right.\n", "\n", "The *shortest path problem*, or how to find the path of minimal distance from a point to another in a network, also belongs in this class. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Fist, let's load libraries we'll need. We'll need `os` and `pandas` for data importation and manipulations, `numpy` and part of `scipy` for computations, `gurobipy` for linear programming using Gurobi, `osmnx` to import map data from [OpenStreetMap](https://www.openstreetmap.org/), `networkx` for operations on networks, parts of `matplotlib` for displays, \n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "ename": "ModuleNotFoundError", "evalue": "No module named 'osmnx'", "output_type": "error", "traceback": [ "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[1;31mModuleNotFoundError\u001b[0m Traceback (most recent call last)", "\u001b[1;32m\u001b[0m in \u001b[0;36m\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[1;32mimport\u001b[0m \u001b[0mosmnx\u001b[0m \u001b[1;32mas\u001b[0m \u001b[0mox\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m", "\u001b[1;31mModuleNotFoundError\u001b[0m: No module named 'osmnx'" ] } ], "source": [ "import osmnx as ox" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "import os\n", "import pandas as pd\n", "import numpy as np\n", "import scipy.sparse as sp\n", "\n", "import networkx as nx\n", "%matplotlib inline\n", "import matplotlib as mpl\n", "import matplotlib.pyplot as plt\n", "from matplotlib.animation import FuncAnimation\n", "import gurobipy as grb\n", "\n", "thepath = 'https://raw.githubusercontent.com/math-econ-code/mec_optim_2021-01/master/data_mec_optim/'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "# A look at our data\n", "\n", "We will be looking at 3 types of data, which we each briefly describe.\n", "\n", "\n", "\n", "## Soviet Planning Problem\n", "\n", "Tolstoi's original data, collected by Schrijver (2002). There are 68 Soviet cities and 10 factories, and 155 links connecting a factory to a city. This yields a sparse $68\\times10$ matrix. Two vectors listing the demand of each city and the supply of each factory are also specified. This is stored in a $69\\times11$ matrix, where we have appended the demand/supply vectors to the right and to the botton of the distance matrix.\n", "\n", "We read the data and display it." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
XArkhangelskYaroslavlMuromBalakhonikhaDzerzhinskKishertSverdlovskArtemovskIledzhkDekonskayademand:
0AgryzNaNNaNNaN709.01064.0693.0NaNNaNNaNNaN2.0
1AleksandrovNaNNaNNaNNaN397.0NaNNaN1180.0NaNNaN4.0
2AlmaznayaNaNNaNNaNNaNNaNNaNNaN81.0NaN65.01.5
3AlchevskayaNaNNaNNaNNaNNaNNaNNaN106.0NaN114.04.0
4BakuNaNNaNNaNNaNNaNNaNNaN1554.0NaN1563.010.0
.......................................
64ShchigryNaNNaNNaNNaNNaNNaNNaN566.0NaN549.04.0
65YudinoNaNNaNNaN403.0757.0999.0NaNNaNNaNNaN0.5
66YamaNaNNaNNaNNaNNaNNaNNaN44.0NaN52.05.0
67YasinovatayaNaNNaNNaNNaNNaNNaNNaN85.0NaN93.06.0
68supply:5.011.58.512.0100.012.015.0314.010.055.0543.0
\n", "

69 rows × 12 columns

\n", "
" ], "text/plain": [ " X Arkhangelsk Yaroslavl Murom Balakhonikha Dzerzhinsk \\\n", "0 Agryz NaN NaN NaN 709.0 1064.0 \n", "1 Aleksandrov NaN NaN NaN NaN 397.0 \n", "2 Almaznaya NaN NaN NaN NaN NaN \n", "3 Alchevskaya NaN NaN NaN NaN NaN \n", "4 Baku NaN NaN NaN NaN NaN \n", ".. ... ... ... ... ... ... \n", "64 Shchigry NaN NaN NaN NaN NaN \n", "65 Yudino NaN NaN NaN 403.0 757.0 \n", "66 Yama NaN NaN NaN NaN NaN \n", "67 Yasinovataya NaN NaN NaN NaN NaN \n", "68 supply: 5.0 11.5 8.5 12.0 100.0 \n", "\n", " Kishert Sverdlovsk Artemovsk Iledzhk Dekonskaya demand: \n", "0 693.0 NaN NaN NaN NaN 2.0 \n", "1 NaN NaN 1180.0 NaN NaN 4.0 \n", "2 NaN NaN 81.0 NaN 65.0 1.5 \n", "3 NaN NaN 106.0 NaN 114.0 4.0 \n", "4 NaN NaN 1554.0 NaN 1563.0 10.0 \n", ".. ... ... ... ... ... ... \n", "64 NaN NaN 566.0 NaN 549.0 4.0 \n", "65 999.0 NaN NaN NaN NaN 0.5 \n", "66 NaN NaN 44.0 NaN 52.0 5.0 \n", "67 NaN NaN 85.0 NaN 93.0 6.0 \n", "68 12.0 15.0 314.0 10.0 55.0 543.0 \n", "\n", "[69 rows x 12 columns]" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "soviet_df = pd.read_csv(thepath+'networks_sovietplanning/distances.csv', sep=',',encoding = \"ISO-8859-1\")\n", "\n", "soviet_df" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plants and cities are the *nodes* of our network. \n", "There are *arcs* directly connecting when there is a railway from the plant to the city. In that case, the entry of the corresponding cell is equal to the length of the arc. When there is none, the corresponding cell shows `NaN`.\n", "We will index the nodes from to `nbPlants+nbCities`, with the plants being indexed by 0,1,...`nbPlants-1`, and the cities being indexed by `nbPlants`, `nbPlants+1`, ... , `nbPlants+nbCities`.\n", "Let us collect the arcs and their corresponding distances.\n", "\n" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "nbPlants = soviet_df.shape[1]-2\n", "nbCities = soviet_df.shape[0]-1\n", "arcsList = []\n", "costslist = []\n", "nodesList = [soviet_df.columns.to_numpy()[plant+1]+' Plant' for plant in range(nbPlants)] + [soviet_df.iloc[city,0]+' city' for city in range(nbCities)]\n", "for plant in range(nbPlants):\n", " for city in range(nbCities):\n", " if not pd.isna(soviet_df.iloc[city,plant+1]):\n", " #arcsList.append((plant,nbPlants+city))\n", " arcsList.append((nodesList[plant],nodesList[nbPlants+city]))\n", " costslist.append(soviet_df.iloc[city,plant+1])\n", "nodesList = [soviet_df.columns.values[plant+1]+' Plant' for plant in range(nbPlants)] + [soviet_df.iloc[city,0]+' city' for city in range(nbCities)]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will create a class `network_problem` to collect and story the relevant information we need: nodes, arcs, and arcs length. We will add methods later." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Number of nodes=78; number of arcs=155.\n" ] } ], "source": [ "class network_problem:\n", " def __init__(self,nodesList,arcsList,costsList):\n", " self.nbNodes = len(nodesList)\n", " self.nbArcs = len(arcsList)\n", " self.nodesList = nodesList\n", " self.arcsList = arcsList\n", " self.costsList = costsList\n", " self.nodesDict = {node:node_ind for (node,node_ind) in zip(self.nodesList,range(self.nbNodes))}\n", " self.arcsDict = {arc:arc_ind for (arc,arc_ind) in zip(self.arcsList,range(self.nbArcs))}\n", " print('Number of nodes='+str(self.nbNodes)+'; number of arcs='+str(self.nbArcs)+'.')\n", "\n", "soviet_pb = network_problem(nodesList,arcsList,costslist) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## NYC subway network\n", "\n", "NYC subway network data: the network data are available on http://www.mta.info. This data is made of two files. The file `nodes.csv` lists the stations: each station is indexed by the line number; each line has the name of the station, and its spatial coordinates. The file `arcs.csv` lists the links between stations: each link specifies the index of the origin, the index of the destination, and the length of the segment. \n", "\n", "We shall store this information in a `network_problem` object called `subway_pb`." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Number of nodes=501; number of arcs=1290.\n" ] } ], "source": [ "arcs = pd.read_csv(thepath+'networks_subway/NYC/arcs.csv', sep=',')#.sort_values(by=['route_id'])\n", "nodes = pd.read_csv(thepath+ 'networks_subway/NYC/nodes.csv', sep=',')\n", "stopnames = [i+' '+j for i,j in zip(nodes.stop_name,nodes.route_id)]\n", "arcsList = [(stopnames[o-1],stopnames[d-1]) for (o,d) in list(zip(arcs.from_stop_nb,arcs.to_stop_nb))]\n", "subway_pb = network_problem(stopnames,arcsList,list(arcs.dis_line)) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Stations' caracteristics are contained in `nodes` dataframe. The relevant information for us will be the name of the station, the horizontal spatial coordinate, and the vertical coordinates." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
stop_nameflowstop_latstop_lonstop_idstop_nbroute_id
0Van Cortlandt Park - 242 St100040.889248-73.8985831011(1)
1238 St100040.884667-73.9008701032(1)
2231 St100040.878856-73.9048341043(1)
\n", "
" ], "text/plain": [ " stop_name flow stop_lat stop_lon stop_id stop_nb \\\n", "0 Van Cortlandt Park - 242 St 1000 40.889248 -73.898583 101 1 \n", "1 238 St 1000 40.884667 -73.900870 103 2 \n", "2 231 St 1000 40.878856 -73.904834 104 3 \n", "\n", " route_id \n", "0 (1) \n", "1 (1) \n", "2 (1) " ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "nodes.head(3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Routes between stations are contained in `arcs` dataframe--the first column is the origin node, the second one is the destination node, and the third one is the length of the arc in meters. " ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
from_stop_nbto_stop_nbdis_linemin_time_elapsedfrom_stop_idto_stop_idroute_id
01851015000.00.18A09112NaN
11982315000.00.18A24125NaN
24482515000.00.18R16127NaN
\n", "
" ], "text/plain": [ " from_stop_nb to_stop_nb dis_line min_time_elapsed from_stop_id \\\n", "0 185 10 15000.0 0.18 A09 \n", "1 198 23 15000.0 0.18 A24 \n", "2 448 25 15000.0 0.18 R16 \n", "\n", " to_stop_id route_id \n", "0 112 NaN \n", "1 125 NaN \n", "2 127 NaN " ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "arcs.head(3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Paris streets\n", "\n", "Next, we gather data on the network of the Paris streets. To to this, we use the `osmnx` package to retrieve data from [OpenStreetMap](https://www.openstreetmap.org/). Doing this is as simple as specifying an address, and calling the `graph_from_address` function, specifying the desired distance (in meters) around this address. " ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "ename": "NameError", "evalue": "name 'ox' is not defined", "output_type": "error", "traceback": [ "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[1;31mNameError\u001b[0m Traceback (most recent call last)", "\u001b[1;32m\u001b[0m in \u001b[0;36m\u001b[1;34m\u001b[0m\n\u001b[0;32m 1\u001b[0m \u001b[0mnyu_address\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m'57 boulevard Saint-Germain, 75005 Paris, France'\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 2\u001b[1;33m \u001b[0mstreets_network\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mox\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgraph\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgraph_from_address\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnyu_address\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdist\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m2000\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mnetwork_type\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;34m'walk'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 3\u001b[0m \u001b[0mox\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mplot_graph\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mox\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mproject_graph\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mstreets_network\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", "\u001b[1;31mNameError\u001b[0m: name 'ox' is not defined" ] } ], "source": [ "nyu_address='57 boulevard Saint-Germain, 75005 Paris, France'\n", "streets_networkn = ox.graph.graph_from_address(nyu_address, dist=2000,network_type = 'walk')\n", "ox.plot_graph(ox.project_graph(streets_networkn))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This network is, however, too big for our trial version of Gurobi. We will therefore Zoom in on a smaller area: the small island at the centre of Paris called Île de la Cité." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "address='Île de la cité, Paris, France'\n", "streets_network = ox.graph.graph_from_address(address, dist=600,network_type = 'walk')\n", "ox.plot_graph(ox.project_graph(streets_network))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As we see, the essential information we can extract from the graph is the list of nodes (`nodes_roads` below) and the list of arcs (`arcs_roads` below)." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "ename": "NameError", "evalue": "name 'streets_network' is not defined", "output_type": "error", "traceback": [ "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[1;31mNameError\u001b[0m Traceback (most recent call last)", "\u001b[1;32m\u001b[0m in \u001b[0;36m\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0marcslist\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m[\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mk\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mlist\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mstreets_network\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0medges\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 2\u001b[0m \u001b[0mcostslist\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m[\u001b[0m\u001b[0mstreets_network\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m\"length\"\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mi\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mj\u001b[0m \u001b[1;32min\u001b[0m \u001b[0marcslist\u001b[0m \u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 3\u001b[0m \u001b[0mstreets_pb\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnetwork_problem\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mlist\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mstreets_network\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mnodes\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0marcslist\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mcostslist\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", "\u001b[1;31mNameError\u001b[0m: name 'streets_network' is not defined" ] } ], "source": [ "arcslist = [(i,j) for (i,j,k) in list(streets_network.edges)]\n", "costslist = [streets_network[i][j][0][\"length\"] for i,j in arcslist ]\n", "streets_pb = network_problem(list(streets_network.nodes),arcslist,costslist) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that we have the structure of the data in mind, let's develop the conceptual tools we shall need in our analysis." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Topology of networks\n", "\n", "### Directed graph\n", "\n", "We start by defining the directed graph on which transportation takes place.\n", "\n", "---\n", "**Definition**\n", "\n", "A (directed) *graph* $\\left(\\mathcal{Z},\\mathcal{A}\\right)$ is a set of *nodes* (cities) $\\mathcal{Z}$, along with a set of *arcs* $\\mathcal{A}\\subseteq\\mathcal{Z}^{2}$ which are pairs $\\left(x,y\\right)$ where $x\\neq y\\in\\mathcal{Z}$.\n", "\n", "---\n", "Our definition does not allow for an arc to have the same origin and destination. Note that for a dense network, $\\left\\vert \\mathcal{A} \\right\\vert = \\left\\vert \\mathcal{Z}\\right\\vert\\left( \\left\\vert \\mathcal{Z}\\right\\vert - 1 \\right)$. For a line, $\\left\\vert \\mathcal{A}\\right\\vert =\\left\\vert \\mathcal{Z}\\right\\vert -1$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Stylized example**. Take the following example, which we shall use for exposition purposes:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Number of nodes=4; number of arcs=5.\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "C:\\Users\\antoi\\Anaconda3\\lib\\site-packages\\networkx\\drawing\\nx_pylab.py:579: MatplotlibDeprecationWarning: \n", "The iterable function was deprecated in Matplotlib 3.1 and will be removed in 3.3. Use np.iterable instead.\n", " if not cb.iterable(width):\n", "C:\\Users\\antoi\\Anaconda3\\lib\\site-packages\\networkx\\drawing\\nx_pylab.py:676: MatplotlibDeprecationWarning: \n", "The iterable function was deprecated in Matplotlib 3.1 and will be removed in 3.3. Use np.iterable instead.\n", " if cb.iterable(node_size): # many node sizes\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADnCAYAAAC9roUQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXxMd/v/8ddkEqIi9rW6SogtfFU1VJuVkNjFTlrctbWC0t5d7tRee6m9i6rWVnvtIWJfG1uKrFUlbag9IrLP74/zQ5UQzMxnluv5eHggM3POe07OXHPNZ875HJ3BYDAghBDCLBxUBxBCCHsiRVcIIcxIiq4QQpiRFF0hhDAjKbpCCGFGUnSFEMKMpOgKIYQZSdEVQggzkqIrhBBmJEVXCCHMSIquEEKYkRRdIYQwIym6QghhRlJ0hRDCjKToCiGEGTmqDiDM71JaJisOJxN3PpXUjBxcnR3xqOBKh1cqU9qlsOp4Qtg0nUxibj+On7vGrB1J7Ey4CEBmTt6d25wdHTAAPtXKMsDbjTrPlVCUUgjbJkXXTiw8cIaxG+PIyMnlYb9xnQ6cHfV8GuRBd68XzZZPCHshwwt2QCu4sdzKznvkfQ0GuJWdy9iNsQBSeIUwMvkizcYdP3eNsRvjClRw/+lWdh5jN8YRk3zNRMmEsE9SdG3crB1JZOTkPtFjM3Jymb0jyciJhLBvUnRt2KW0THYmXHzoGO7DGAywPf4il9MyjRtMCDsmRdeGrTicnO9tf4xvQfbVv+78/9L6qVzd9eN999MBK47kvxwhxOORomvD4s6n3nNY2JPIyMkjLuWGkRIJIaTo2rDUjBwjLSfbKMsRQkjRtWmuzsY5ItDV2ckoyxFCSNG1aR4VXCns+OBfsc6pMIbsu1+Q5d68+sD7OTs64FGxmEnyCWGPpOjasJBXKud7W6FyL3Pz1E4MebncOn2YzHMnHng/AxBSL//lCCEejxRdG1bGpTDeVcui091/W8mAPtxKOsS5aZ25eXIHRdy97ruPTge+1crKJDhCGJHMvWDjjp+7RudvDnAr+/FPkCji5MBPfRriWVkmvxHCWKTTtXF1nivBp0EeFHF6vF+1Ljebyhd/oaaM5wphVFJ07UB3rxf5NKg6RZz0DxxquIchjyJOesJb1ICk3fTo0YOcHOMceiaEkKJrN7p7vUhrl9NUK5pJYUcHnP91VIOzowOFHR3Qp5xkYI1cer1ZlXXr1nHt2jU6d+5MdrYcqyuEMciYrh3IyMhgwIABzJ8/n06dOjHr2wWsOJJMXMoNUjOycXV2wqNiMULqVWZHxHo+//xzoqOj0el0ZGZm0rFjRwCWLVtG4cLypZoQT0OKro2LjY2lZcuWnD17luzsbMLDwxk1alS+98/Ly6N+/fqEh4fTtm1bALKysujatSvp6emsXLmSIkWKmCu+EDZHhhdsWE5ODl5eXpw+ffrO8ICj48PPUnNwcGDUqFF89tln5OVp8zYUKlSIpUuXUqJECVq1akV6errJswthq6To2jBHR0d27dpF9erVcXBwQKfTodfrH/m44OBgihYtyrJly+5Z1o8//sizzz5LUFAQaWlppowuhM2SomvjqlevTnp6OnPmzOHNN9+kSpUqj3yMTqdj9OjRDB8+/J4jF/R6Pd999x3u7u4EBgaSmppqyuhC2CQZ07Vxc+fOZc2aNWzevPmxHmcwGPD29qZ379689dZb99yWl5fHwIEDiY6OZvPmzZQsWdKYkYWwaVJ0bVhGRgbu7u6sXLmSBg0aPPbjd+3axdtvv018fDxOTvfONGYwGBg6dCg7duxg69atlC5d2lixhbBpMrxgw7766ivq1av3RAUXuDMcMX/+/Ptu0+l0TJkyhcDAQHx9ffn777+fNq4QdkE6XRuVnp5OlSpV2Lx5M3Xq1Hni5Rw4cICOHTuSkJCAs7PzfbcbDAZGjBjB8uXL2bZtGxUrVnya2ELYPOl0bdSsWbNo3LjxUxVcAC8vLzw9Pfnmm28eeLtOp2PkyJF069YNb29vkpPlempCPIx0ujboxo0buLm5sX37dmrUqPHUyzty5AgtWrQgKSmJZ555Jt/7TZ48mblz57Jt2zZeeOGFp16vELZIOl0b9OWXX9KkSROjFFyAevXq0bBhQ+bMmfPQ+w0bNoywsDB8fHw4ffq0UdYthK2RTtfGXL16FXd3d/bv34+7u7vRlnvixAn8/f1JSkqiWLGHT/f41VdfMXbsWCIjI6latarRMghhC6TTtTFffPEFrVu3NmrBBahVqxb+/v7MmDHjkfft27cvI0aMwM/Pj9jYWKPmEMLaSadrQy5dukS1atU4fPgwL774otGXHx8fT+PGjUlMTKREiUdfTWLhwoV8+OGHREREULt2baPnEcIaSadrQyZOnEinTp1MUnABqlWrRnBwMFOnTi3Q/bt3787UqVNp0qQJR48eNUkmIayNdLo24vz589SsWZOYmBieffZZk63n9OnTNGjQgPj4+AKfhbZ69Wr69evHunXrnvhEDSFshRRdGzFo0CB0Oh3Tpk0z+br69u1LyZIlGT9+fIEfs379enr16sWaNWto1KiRCdMJYdmk6NqA5ORk6tSpw8mTJ6lQoYLJ13fu3Dnq1KlDbGws5cuXL/DjIiIi6NGjBytWrODNN980YUIhLJcUXRvQv39/XF1dmTBhgtnWGRYWhl6vL/D47m1RUVF07tyZJUuW4O/vb6J0QlguKbpW7vfff6d+/fokJCSYdaavlJQUatasya+//vrYY8i7d++mffv2/PDDDzRr1sxECYWwTFJ0rVyvXr2oXLnyQ697ZioffPAB6enpzJo167Efu3//ftq0acO3335Ly5YtTZBOCMskRdeKJSQk8Prrrxf4uFlju3jxIh4eHhw5cuSJ5lqIjo4mODiY2bNn0759exMkFMLyyHG6VmzkyJEMHjxYScEFKFu2LP369WP06NFP9Pj69esTERHBe++9x9KlS42cTgjLJJ2ulTp58iR+fn4FmgvBlG7P9XDgwAHc3NyeaBknTpygadOmjB8/ntDQUCMnFMKySKdrpYYPH86wYcOUFlyAkiVLEhYWxsiRI594GbVq1SIqKopPPvmEefPmGTGdEJZHOl0rdOzYMYKCgh45v625pKam4ubmxo4dO55qOsnExEQCAgL473//y4ABA4yYUAjLIUXXCrVq1YqAgADCwsJUR7ljwoQJHD58mGXLlj3Vcn7//Xf8/PwYNGgQgwcPNlI6ISyHFF0rc/DgQUJCQkhMTHzgNctUuXnzJlWqVCEiIuKpLxF09uxZ/Pz86NOnDx9++KGREgphGaToWpnAwEDatWtH3759VUe5z7Rp09ixYwdr1qx56mX9+eef+Pv7061bN8LDw42QTgjLIEXXiuzevZvQ0FDi4+MpVKiQ6jj3ycjIwM3NjdWrV/Pqq68+9fLOnz+Pv78/7dq1Y9SoUeh0OiOkFEItKbpWwmAw4Ovry9tvv83bb7+tOk6+5syZw9q1a9m0aZNRlnfx4kUCAgIIDAxkwoQJUniF1ZNDxqxEVFQUKSkpdO/eXXWUh+rduzexsbHs3bvXKMsrW7YsUVFRbNu2jSFDhiA9grB20ulaAYPBQKNGjQgLC6NLly6q4zzSvHnzWLRoEVFRUUZb5rVr12jWrBn16tVj5syZODhIvyCsk+y5VmDTpk3cuHGDTp06qY5SIKGhoZw7d86oRbdEiRJs2bKFX3/9lT59+pCbm2u0ZQthTlJ0LZzBYCA8PJyRI0daTXfn5OTE8OHDCQ8PN+pwgKurK5s2beK3336jZ8+e5OTkGG3ZQpiLdbyK7diaNWswGAy0bdtWdZTH0qVLF65evUpERIRRl+vi4sKGDRs4f/483bt3Jzs726jLF8LUZEzXguXl5VGnTh3GjRtHixYtVMd5bMuXL2fixIkcOnTI6EcdZGRk0L59e5ydnVmyZIlFHkInxINIp2vBli1bRtGiRQkODlYd5Ym0b9+erKws1q5da/RlOzs7s2rVKnJzcwkJCSEzM9Po6xDCFKTTtVA5OTnUqlWLGTNm0KRJE9VxntjatWsJDw/n6NGjJhmTzs7Oplu3bqSmprJ69WqKFCli9HUIYUzS6VqoxYsXU65cOQICAlRHeSotW7akcOHCrFixwiTLd3JyYvHixZQuXZoWLVpw8+ZNk6xHCGORTtcCZWdn4+Hhwfz5823iUuUREREMHjyYEydOoNfrTbKO3Nxc/vOf//Dbb7+xYcMG5fMMC5Ef6XQt0Pfff8/LL79sEwUXoGnTppQuXZrFixebbB16vZ558+ZRvXp1AgMDuX79usnWJcTTkE7XwmRmZuLu7s6yZcvw8vJSHcdotm/fzjvvvENsbCxOTk4mW4/BYCAsLIyDBw8SERFByZIlTbYuIZ6EdLoW5ptvvsHT09OmCi6Ar68vL7zwAj/88INJ16PT6Zg+fTpvvPEGfn5+XLp0yaTrE+JxSadrQdLT03F3d2fdunXUq1dPdRyj27dvH126dCEhIYHChQubdF0Gg4FPP/2UdevWERkZSfny5U26PiEKSjpdCzJnzhy8vLxssuACNGrUiBo1apjl4pM6nY6xY8cSEhKCj48PKSkpJl+nEAUhna6FSEtLw83NjcjISGrVqqU6jslER0fTunVrkpKSzHZM7bhx45g/fz5RUVFUrlzZLOsUIj/S6VqIGTNm4Ovra9MFF6B+/fq8+uqrzJ0712zr/Pjjj+nXrx/e3t6cOXPGbOsV4kGk07UA169fx83NjT179lCtWjXVcUwuJiaGpk2bkpSUhIuLi9nWO3PmTCZPnsy2bduoUqWK2dYrxD9Jp2sBpk6dSnBwsF0UXABPT0+8vb2ZOXOmWdf73nvv8cknn+Dj40N8fLxZ1y3EbdLpKnb58mWqVq3KL7/8wssvv6w6jtnExsbi7e1NYmIixYsXN+u658+fz//+9z+2bt1KjRo1zLpuIaTTVWzy5MmEhITYVcEFqF69Os2aNWPatGlmX3fPnj2ZOHEiAQEBxMTEmH39wr5Jp6vQ33//TfXq1Tl27BjPPfec6jhml5SUhJeXFwkJCZQqVcrs61++fDkDBw5k48aNNnuYnrA80ukqNH78eLp162aXBRfAzc2Ntm3bMmXKFCXr79ChA3PnzqV58+YcPHhQSQZhf6TTVeTPP//E09OTEydOULFiRdVxlPnjjz+oV68ecXFxlC1bVkmGDRs20LNnT1atWkXjxo2VZBD2Q4quIu+++y5FihRh8uTJqqMoZwnbYsuWLXTv3p1ly5bh4+OjLIewfVJ0FbCE7s6S/PXXX9SqVYuTJ08q7fq3b99Op06dWLx4sdVPHi8slxRdBd555x3KlSvH2LFjVUexGO+//z7Z2dnMmDFDaY7du3fTvn17vv/+e4KCgpRmEbZJiq6Z3f7GPjExUeZ6/YfbR3IcPXqU559/XmmWAwcO0KpVK7755htat26tNIuwPXL0gpmNGjWKsLAwKbj/Uq5cOfr06cOYMWNUR8HLy4tNmzbRt29fk13bTdgv6XTN6PZZWElJSbi6uqqOY3Es7ey848eP06xZM6ZMmULXrl1VxxE2QjpdMxoxYgRDhw6VgpuP0qVLM3DgQEaNGqU6CgB16tQhMjKSDz74gAULFqiOI2yEdLpmEhMTQ2BgIElJSRQtWlR1HIt1e8a13bt34+HhoToOAHFxcTRp0oTPPvuMd955R3UcYeWk6JpJmzZt8Pb2ZsiQIaqjWLxx48YRExPDkiVLVEe5IykpCX9/fz788EPeffdd1XGEFZOiawbR0dG0adOGxMREs10twZrdvorG1q1bqV27tuo4d5w5cwY/Pz8GDhwob57iiUnRNYOgoCBatGjBgAEDVEexGlOmTGHv3r2sWrVKdZR7nDt3Dj8/P3r37s1HH32kOo6wQlJ0TWzfvn107dqV+Ph4k18B15bcvjLy2rVreeWVV1THucdff/2Fn58fXbt2JTw8HJ1OpzqSsCJSdE3M39+frl270rt3b9VRrM7MmTPZtGkTGzZsUB3lPhcuXMDf35/WrVszZswYKbyiwOSQMRPavn07f/zxB6GhoaqjWKV33nmHEydOsH//ftVR7lO+fHm2b9/Ohg0b+OCDD5DeRRSUFF0TMRgMpKSkEBMTg5OTk+o4Vqlw4cKcPHnSYq/gW7ZsWaKioti5cyeDBg2SwisKRIYXTCQnJweDwSAF1wiysrLQ6/Xo9XrVUR4oLS2NjIwMSpQogV6vl6EG8VBSdIUQwoxkeEEIIcxIiq4QQpiRo+oANun4cVi8GG7dgs6doVEj1YmsX24u6HTgYCV9Ql6e9re15BVmI3uEsS1fDqGhULgwNGgA3brBjRuqU1mvrCztb73eugqYg4N15RVmI1+kGdvNm/DPWcQ8PWHhQu1v8Xji4uC77yAyEry9oVcvqF1b63odHLTO1xKdOAFr1oCbG7z6KlSpAgbD3e5Xp7v7B+DKFZg3DwoVgpo1oXFjcHZWl1+YlLwVP63c3LsvJrhbcKOitEL7/PN3uzVRcBcuwMSJ2rbduBEqV4bt27Xb9Pq7Bevf2x/g+nXzZv2no0ehXz/IzobNm2HaNO3nOp2W+3bH/s83jAsXtL9Pn4bx42H6dO15CZsknW5B3O5Sbh8n+qBO6+pVSEmBGjW0/+/YAb//DqVLwxdfaC8k6XYLbskSSEyE//wHKlXSfhYbq23LQoXgrbe04Zt/u3wZgoK0N7oRI8Cc1zjLyoLPP4fMTBg3TvtZWhp8/TUsWKB1vG5uWrdeuzZUr64NQ2Vlac8JtO8BXnwR/vhDul0bJZ1uQdzuUkAbPrjdaf39N8yeDc2bQ0AA7Nt39zE+PtCzJ7RqBdWqwYEDSqJbrawsrfBUqqR1jaBtx4kToV49GDVKK7oDB8LevXcfV6iQVpirV7/7+zBXX5GWpu0fgYHa/7OztU8+778Pe/bAyJFQt642/BAeDlOm3M18u1uPj9eKc06OeTILs5OjFwri4kWYNAkOHdLGFlu2hGHDtK7Ezw8mTNC62Ph47aNt8eL3Pj4z8+4L32Cw3LFIS1Kq1N1hAicnrQhPnaoV0h494KeftE8c69fDjz9q279YMe3Ps8+CoyM0aaI93lzbvEQJrTs/f/5u7t9/194gkpKgTh2tw23TBj78UPsUBFqBdXTUHjt7tvam4uIi+4qNkuGFghg7Fi5dgt69teGDuDjtBeXuro3hLVmijTeWKaN1WS+9pL2Ali3TOpxnnoG5c+G551Q/E+tx48bdDvHNN7XCVb48bN0Kq1drb4CVK0NICLRvD7cnhzcYtC/epk+HlSvvfmw3lx9+gHXrtC/DqlSBpk21DPHxcOoUxMRowyTnz8Po0fDGG9rjtm3T9p369WH4cG34SoquTZKi+yg5OdqLfsoUaNhQ+xiYmal92zxrllZ4g4O1F36ZMtpjcnO1L1G2bYOOHcHLS+1zsGZXrmjDB/7+2pvXP+3ZoxXW5s214gbaR/wZM7Sx3cmTtd+XOQ/dyszU3hSio7V9Z+TI+z/53Ha7qP74I0REaIca3n4ewmZJ0X2Ua9e0TrdJk7sviOxs7QSIqlXhn1f2zcvTXkgWOjGLTXjUNj5/HoYMga5dtWEgcxfdh7mdHe6e6LF5s/bFX2CgdtRD1arg4SEdrg2TMd1HcXXVOtjFi+8W3QMH4PXXtRfN7UN7rO3gfWv1z238zwJ886Z2Ykrx4trvpGXL+++v2oOyvPYarFqljf3+9JN2BMyXX8qRLjZMOt2CuH4d3n5bG5s7dUobUpg4UTv8R1iG5GQIC9O+aLtyRRtbbdtWOyRLCAsiRbegbtyAI0fghRe04yiF5Vq9Gs6ehe7d7x4hIISFkKIrhBBmZEEDXkIIYfuk6D6B3Nxc8v59vr9QypZ+Jzk5OWRkZHBDZqezSVJ0H1NGRgbR0dE4WNK34gK9Xs/BgwfJzMxUHeWpOTo6sm3bNtzc3IiOjlYdRxiZVI7HcOnSJSpXrkzlypVVRxEPULFiRZ577jmuq5xlzEiCg4P5+uuvCQoK4oDM22FT5Iu0xxAWFoajoyNffPGF6igiH/3798fV1ZUJEyaojmIUmzZt4q233mLVqlU0btxYdRxhBFJ0C+jcuXPUrVuX2NhYypUrpzqOyEdycjKenp6cOnWKChUqqI5jFJGRkXTp0oVly5bh6+urOo54SlJ0C6hv376ULFmS8ePHq44iHmHw4MEATLs9gbgN2LFjBx07dmThwoU0lfkZrJoU3QI4ffo0DRo0ID4+ntJysL3FO3/+PDVq1CAmJsamxt/37t1L27ZtmT9/PsHBwarjiCckRbcAevbsyfPPP8/IkSNVRxEF9N///pfU1FTmzJmjOopRHTx4kFatWvHVV1/Rpk0b1XHEE5Ci+wjx8fE0btyYxMRESpQooTqOKKBLly5RrVo1oqOjeemll1THMaojR44QFBTEjBkz6NChg+o44jHJIWOPMHLkSIYMGSIF18qUKVOGAQMGMHr0aNVRjK5evXpEREQQFhbGokWLVMcRj0k63Yc4ceIE/v7+/Pbbb7i4uKiOIx7TtWvXcHd3Z9++fbi7u6uOY3QnT56kadOmjBkzhp49e6qOIwpIiu5DtG/fnkaNGjF06FDVUcQTGjNmDHFxcSxcuFB1FJOIj48nICCA8PBw+vTpozqOKAApuvk4cuQILVq0ICkpiWf+fZkYYTVu3LiBm5sbUVFR1KxZU3Uck/jtt9/w9/dn2LBhvPfee6rjiEeQopuPFi1aEBgYyMCBA1VHEU9p0qRJHDp0iOXLl6uOYjJnzpzBz8+Pd999Vz6ZWTgpug9w4MABOnbsSGJiIoXlygNWLz09nSpVqrBp0ybq1q2rOo7JnDt3Dn9/f95++20++eQT1XFEPqToPkCTJk3o0KGDjJHZkOnTpxMZGcnatWtVRzGplJQU/Pz86NSpE8OHD0cnF7i0OFJ0/2XXrl307NmTuLg4nJycVMcRRpKRkYG7uzsrV66kQYMGquOY1IULFwgICKBly5aMHTtWCq+FkaL7DwaDAW9vb3r37s1bb72lOo4wsq+++opVq1YRERGhOorJXbp0iSZNmuDn58fkyZOl8FoQOTniHyIjI7lw4QLdunVTHUWYQM+ePUlISGDPnj2qo5hcmTJl2LZtG7t27SIsLMxmrqphC6TT/f8MBgMNGzZk8ODBdO7cWXUcYSLz589nwYIFbN++3S66v+vXr9O8eXNq1arF3Llz5YonFkB+A//fxo0buXnzJh07dlQdRZhQjx49+Ouvv4iKilIdxSyKFy9OREQE8fHx9O7dm9zcXNWR7J4UXbQuNzw8nFGjRkknYOMcHR0ZMWIE4eHh2MuHvGLFirFx40bOnj1LaGgoOTk5qiPZNakwwOrVq9HpdDJVnp3o1KkTqampbNq0SXUUsylatCjr16/n8uXLdO3alezsbNWR7Jbdj+nm5uZSp04dJkyYIBND25GVK1cybtw4fvnlF7sY270tIyODDh06oNfr+emnn+TkHwXsvtNdtmwZxYoVIygoSHUUYUZt27YlNzeXn3/+WXUUs3J2dmblypXo9XratWtHRkaG6kh2x6473ZycHGrWrMmsWbMICAhQHUeY2fr16/n44485fvy43Y3lZ2dnExoayuXLl1mzZo1M6mRG9rWn/cvChQupWLEi/v7+qqMIBYKDgylatKhNT4STHycnJxYuXEiFChUIDg4mLS1NdSS7YbedblZWFh4eHixYsIA33nhDdRyhyNatWxk4cCAnTpzA0dFRdRyzy83NpW/fvsTFxbFx40ZcXV1VR7J5dtvpzp8/Hzc3Nym4di4gIIBy5cqxePFi1VGU0Ov1fP3113h6etK0aVOuXbumOpLNs8tO9/bkJytWrOC1115THUcotnPnTnr16mXXkxwZDAaGDBnCnj172LJlC6VKlVIdyWbZZaf7zTffULduXSm4AgBvb29efvllvv/+e9VRlNHpdEydOhV/f398fX25ePGi6kg2y+463fT0dNzc3NiwYQP/93//pzqOsBAHDhygU6dOJCQk2PWxqwaDgc8++4xVq1axbds2KlSooDqSzbG7Tnf27Nk0atRICq64h5eXF7Vr1+bbb79VHUUpnU7H6NGj6dy5Mz4+Pvz555+qI9kcu+p07eEiheLJHTlyhJYtW5KUlESRIkVUx1Fu4sSJfP3110RFRfH888+rjmMz7KrTnT59OgEBAVJwxQPVq1eP1157jTlz5qiOYhE+/PBD3nvvPby9vfn9999Vx7EZdtPpXrt2DXd3d/bu3UvVqlVVxxEW6tdff6VJkyYkJSXh4uKiOo5FmDNnDuPHjycyMhJ3d3fVcaye3XS6X3zxBS1atJCCKx6qdu3a+Pr6MmPGDNVRLEb//v357LPP8PX1JTY2VnUcq2cXne6lS5eoVq0a0dHRvPTSS6rjCAsXHx9P48aNSUpKonjx4qrjWIwffviBjz76iC1btlCrVi3VcayWXXS6kyZNokOHDlJwRYFUq1aN4OBgpk6dqjqKRQkNDeWLL76gSZMmHDt2THUcq2Xzne6FCxeoUaMGx48fp3LlyqrjCCtx+vRpGjRoQEJCgpyd9S8rV65kwIABbNiwgfr166uOY3VsvugOGTKEvLw8vvzyS9VRhJXp27cvpUuX5vPPP1cdxeKsW7eO3r178/PPP9OwYUPVcayKTRfd5ORkPD09OXXqlJxZIx7buXPnqFu3LrGxsZQrV051HIuzefNmQkNDWblypUwc9RhsuugOGDAAFxcXJk6cqDqKsFIDBw6kUKFCTJkyRXUUi7Rt2za6dOnC0qVL8fPzUx3HKths0T1z5gyvvPIK8fHxlClTRnUcYaVSUlKoWbMmJ06coFKlSqrjWKSdO3fSoUMHfvzxRwIDA1XHsXg2W3R79+5NpUqVGD16tOoowsoNGzaMjIwMZs6cqTqKxdq3bx9t2rThu+++o0WLFqrjWDSbLLqJiYk0bNiQxMRESpYsqTqOsHIXL17Ew8ODI0eO8MILL6iOY7EOHTpEy5YtmTt3Lm3btlUdx2LZ5HG6I0eOZNCgQVJwhVGULVuWfv36MWbMGNVRLFqDBg3YtGkT/fv356efflIdx2LZXKd76kvwf8kAAAzxSURBVNQpfHx8SEpKkus9CaO5cuUKVatW5eDBg1SpUkV1HIv266+/EhgYyMSJE+nevbvqOBbH5jrd4cOHM2zYMCm4wqhKlSpFWFgYo0aNUh3F4tWuXZvIyEg++ugjvvvuO9VxLI5NdbrHjh2jefPmJCUlUbRoUdVxhI1JTU3Fzc2NXbt24eHhoTqOxUtISCAgIIBPPvmEfv36qY5jMWyq6LZu3Ro/Pz8GDRqkOoqwURMmTODo0aMsXbpUdRSrcPr0afz8/Hj//fcJCwtTHcci2EzRPXToEO3btycxMRFnZ2fVcYSNunnzJm5ubkRERODp6ak6jlX4448/8PPzo3///gwbNkx1HOVspug2a9aMNm3ayMcYYXJTp05l165drF69WnUUq5GcnIy/vz+hoaF8+umnquMoZRNFd8+ePfTo0YP4+HgKFSqkOo6wcbdu3cLd3Z2ff/6ZV155RXUcq5GSkkJAQAAhISGMGDECnU6nOpISNlF0fX196dGjB7169VIdRdiJ2bNns379ejZu3Kg6ilX5+++/CQgIICgoiHHjxtll4bX6Q8aioqJITk4mNDRUdRRhR3r37s2pU6fYv3+/6ihWpVy5cmzfvp0tW7YwdOhQbKDne2xW3ekaDAZef/113n33Xbp166Y6jrAz8+bNY8mSJURGRqqOYnWuXr1Ks2bNePXVV5k+fToODlbf/xWYVT/TzZs3c/36dTp37qw6irBDoaGhnDlzhh07dqiOYnVKlizJ1q1bOXr0KH379iUvL091JLOx2qJrMBgIDw9n5MiR6PV61XGEHXJycmLEiBGEh4fb5cfkp+Xq6kpERASJiYn06tWL3Nxc1ZHMwmqL7s8//0xOTg7t2rVTHUXYsS5dunD58mW2bNmiOopVcnFxYePGjfz555/06NGDnJwc1ZFMzirHdPPy8qhbty5jx46lZcuWquMIO7d8+XImTZrEwYMH7fLbeGO4desW7dq1w8XFhcWLF+Pk5KQ6kslYZae7fPlyihQpIpMlC4vQvn17MjMzWb9+veooVqtIkSKsWbOGrKwsQkJCyMzMVB3JZKyu083NzaVWrVpMmzZNLg0iLMbPP//M8OHDOXLkiF19E29sWVlZdO3alfT0dFauXEmRIkVURzI6q9s7Fi9eTJkyZWjatKnqKELc0apVK5ycnFi1apXqKFatUKFCLF26lBIlStCqVSvS09NVRzI6q+p0s7Oz8fDwYN68efj4+KiOI8Q9Nm/ezNChQ4mJiZEjap5Sbm4uvXr14uzZs6xbtw4XFxfVkYzGqjrdBQsW8NJLL0nBFRYpMDCQkiVLyrSPRqDX65k/fz5ubm40a9aM1NRU1ZGMxmo63czMTKpWrcrSpUtp2LCh6jhCPND27dvp06cPsbGxODo6qo5j9fLy8hg4cCDR0dFs3rzZJq57aDVFd9asWWzYsEEmGBEWz8/Pj44dO5KWloaDgwPvv/++6khWzWAwMHToUHbs2MHWrVspXbq06khPxSqK7q1bt3Bzc2Pt2rUylZ6waJmZmXz00UdMnz4dvV5P/fr12bdvn+pYVs9gMPDxxx+zceNGIiMjKVeunOpIT8wqxnTnzJlDgwYNpOAKi9emTRtmzpxJXl4e2dnZdnNqq6npdDrGjRtH27Zt8fX1JSUlhevXr/P5559b3SnYFj/olJaWxsSJE9m6davqKEI80rRp02jZsiVnz54lMzPTLk5rNRedTsfIkSMpVKgQb775Jk5OTsTFxeHn54eXl5fqeAVm8Z3ujBkz8PHxoXbt2qqjCPFI1apVIyYm5s78zhcuXFCcyPYMHjyYW7duERsbi06nY/78+aojPRaLGNO9lJbJisPJxJ1PJTUjB1dnRzwquBJYtTgNPKuze/duueS1sDqzZs1i8+bNrFu3Lt99vMMrlSntUlh1VKvSqVMnVq5ceWfoplixYly9evWeY6MteXsrLbrHz11j1o4kdiZcBCAz5+6cms6ODmTn5FDy1l/MG9qJOs+VUBVTiCf2qH3cAPhUK8sAbzfZxwvo5MmTfPvttyxdupQrV66QlZXFwoUL6datm1Vsb2VFd+GBM4zdGEdGTi4PS6ADnJ30fBrkQXevF80VT4inVuB9XAfOjrKPPy6DwcDJkycZM2YMISEhZFSubxXbW0nR1XbGWG5lF3y2+CJODnwaVF12SmEVZB83L2va3mb/Iu34uWuM3Rj3WBsH4FZ2HmM3xhGTfM1EyYQwDtnHzeth2zt5di9unTn2wMep2t5mL7qzdiSRkfNkxy5m5OQye0eSkROJB0pNhWefhffeU53E6tjFPn7sGDRsCDVrgqcn/PSTsijWtr3NWnQvpWWyM+FivuMt1w+sIHlmKGe/6MCfX/e97x3KYIDt8Re5nGa7ExxbjPBw8PZWncJi7dmzh6VLl9538kN++/j1/cu5uPrze352ZetXXNn61T0/s5p9/Jln4Icf4ORJ2LwZBg+Ga6brGM+ePcuMGTNIS0u75+cPqymX1k0hN/UiF1eM4uyUEK4fWHHffVRsb7MW3RWHk/O9LftyMjcOr6fCW1N5/v3llO84Csfi5e+7nw5YcST/5Yh/+eEHrROpUwd69CjYYw4fhgsXQOYszteiRYvo3r07L774IosXL75TfPPbx4vW8ObWb4fJy9TmhzXk5ZIet4eiNe5/Y1Oyjz/uflK1Kri7a/+uVAnKlYOLF00W7+DBgwwZMoRKlSrx+eef3ym+D6spZVoORe9alrIhn/H80BUU9wp54P3Mvb3NekZa3PnUew7huIeDA4bcbLIvn0P/THEcS9xfcAEycvL4ce02Tq740oRJrU/RokUZNWrUvdeWOnkSxo6FvXuhTBm4cgUWLYJJk+5fgJsbrFgBeXkwdCj8+CNs2/bAdeXk5PDRRx+RlZVlomdj+fbs2UNubi7Jycm89dZb9O3bl0WLFhF367kH7uOOxctRqEIV0hP241Lbn4w/YtA5Fabws/cff27qffx///vfvXMXPMl+8k+HDkFWFlSpcs+Pc3JyWL9+PVFRUU+dOTExEQcHB27cuMHw4cMZMWIEQ4cO5VbdjvnXlALKyMkjLuXGU2csKLMW3dSM/E+JdCpZiVL+73B9z2IuXfwD55frUdLvPzgWu39GoUIuJXAr5WbKqFbH2dn5/osiRkVBSIj2QgIoVQq6ddP+5Gf2bAgKgueey/cuOp2OKlWqkJ2dbYTk1unw4cN3/q3T6ShfvjwVKlQgNTb/fbxoDW9uxu7CpbY/N0/tfGCXe5sp9/H7Lvr4JPvJbSkpWme8YAH86zJFOp2OMmXK4Ob29M/jxo0bdy6D5ODggIuLCy+//DL7H1JTHkdqhvn2ZbMWXVfnh6+uaE0fitb0IS8zncubZ3Jtx/eUaTn0vvvVqV6VsE4dTRXTdhgM2kGJ//SoDmb/fti9Wyu+aWlaB+PiAuPH37mrXq+nf//+Jg5v2WJjY9m/fz9vvPEGkyZNokGDBgAs/v1ovo95xqMxV6PmkZN6ifSE/VTsMTnf+5p1H3+S/QS0L1uDg2HMGHjA3Ad6vZ7GjRvTuHHjp464fPnyOxcxmDhxIu3atcPBwYGTP+W/vYH7n1c+XJ3Nd/VhsxZdjwquFHY8/8CPA9mXk8m5cRnnyjXQOTqhcyzEg0bHnR0d8KhYzBxxrZ+/P7RtC0OGQOnS2sfGR3Uwixbd/ff330N09D0FV2gGDRpEr169ePXVV+/5+cP2cf0zxSn8fG0ub5yGY/HyOJV58KcJs+/jT7KfZGVpjwkNhQ4dzBDRnw0bNhAYGHjPhT8ftr0B9EVLkHPt/EOXbe7tbdYv0kJeqZzvbYbcbK7t/J5z07uSPKMHeenXKeEdev/9gJB6+S9H/EPNmvDpp9pRCHXqgEymbTQeHh73FVx4+D4O2hBDxpljFK2Z/9CC2ffxJ9lPli2DXbu0N+a6dbU/xx58PKwxlCpViubNm993peVHbe/iXh24vu8nzk7txPWDD75oqLm3t9nPSOvzYzRbYy889DS9/Oh0EFijPHO71zd+MCGMRPZx87K27W32kyPe9XHD2fHJrpTq7KhngI98gSYsm+zj5mVt29vsRbfOcyX4NMiDIk6Pt2rtPGkPPCvLTEzCssk+bl7Wtr31I0aMGGHWNQKelUtQoogT+09fIfcRnwl0OijipJeJQIRVkX3cvKxpeyudTzcm+RqzdySxPf4iOrSDlG+7Pfelb7WyDPBxk3d/YZVkHzcva9jeFnHliMtpmaw4kkxcyg1SM7JxdXbCo2IxQuqpn+VdCGOQfdy8LHl7W0TRFUIIe2HxF6YUQghbIkVXCCHMSIquEEKYkRRdIYQwIym6QghhRlJ0hRDCjKToCiGEGUnRFUIIM5KiK4QQZiRFVwghzEiKrhBCmJEUXSGEMCMpukIIYUZSdIUQwoyk6AohhBn9P+pRLgyH5UN+AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "example_graph = nx.DiGraph()\n", "arcslist = [('s','u'),('u','v'),('s','v'),('u','t'),('v','t')]\n", "costslist = [3,1,4,2,2]\n", "example_graph.add_edges_from(arcslist)\n", "\n", "example_pb = network_problem(list(example_graph.nodes),arcslist,costslist) \n", "\n", "pos = {'s':np.array([0,0]),'u':np.array([1,0.5]),'v':np.array([2,0]),'t':np.array([3,0])} \n", "nx.draw_networkx_edge_labels(example_graph,pos,edge_labels={e:('c='+str(l)) for (e,l) in zip(arcslist,costslist)},font_color='red')\n", "nx.draw(example_graph,pos,with_labels=True)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Gradient matrix\n", "\n", "Next, we define the gradient matrix, which encodes the information contained in the directed graph.\n", "\n", "---\n", "**Definition**\n", "\n", "\n", "We define the *gradient matrix* (also called an `edge-node matrix') as the matrix with general term $\\nabla_{ax}$, $a\\in\\mathcal{A}$, $x\\in \\mathcal{Z}$, such that\n", "\n", "\\begin{align*}\n", "\\nabla_{ax} & =-1\\text{ if }a\\text{ is out of }x\\text{,}\\\\\n", "& =+1\\text{ if }a\\text{ is into }x\\text{,}\\\\\n", "& =0\\text{ else.}\n", "\\end{align*}\n", "\n", "---\n", "\n", "Hence, if $f\\in\\mathbb{R}^{\\mathcal{Z}}$, then $\\nabla f\\in\\mathbb{R}^{\\mathcal{Z}}$, and $\\left( \\nabla f\\right) _{xy}=f_{y}-f_{x}$.\n", "\n", "We shall denote $\\nabla^{\\intercal}$ the transpose of the gradient matrix. It is the network analog of the $-\\operatorname{div}$ differential operator.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "From the list of arcs, now we can easily construct our $\\nabla$ matrix, as the first component of each element of `arcs` is the origin node (which will have $-1$ entry in $\\nabla$) and the second component is the destination nodes (which will have $+1$ entry). \n", "An important feature of the matrix $\\nabla$ is *sparsity*: for each row (arcs), there are only two nonzero terms. We should recognize this sparse structure and code $\\nabla$ as a sparse matrix. In `scipy`, this is done using `csr_matrix`. \n", "We implement this as follows:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "def Nabla(self):\n", " data = np.concatenate([-np.ones(self.nbArcs),np.ones(self.nbArcs)])\n", " arcsIndices = list(range(self.nbArcs))\n", " arcsOrigins = [self.nodesDict[o] for o,d in self.arcsList]\n", " arcsDestinations = [self.nodesDict[d] for o,d in self.arcsList]\n", " theNabla = sp.csr_matrix((data, (arcsIndices+arcsIndices, arcsOrigins+arcsDestinations)), shape = (self.nbArcs,self.nbNodes))\n", " return(theNabla)\n", "\n", "network_problem.Nabla = Nabla" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we want to visualize, we need to convert it to a `np.array` matrix using the `.toarray()` method. \n", "\n", "**Stylized example**. In the case of our example, the matrix looks like the following:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[-1., 1., 0., 0.],\n", " [ 0., -1., 1., 0.],\n", " [-1., 0., 1., 0.],\n", " [ 0., -1., 0., 1.],\n", " [ 0., 0., -1., 1.]])" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "example_pb.Nabla().toarray()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Take the first line in the example above, for instance, which is associated with the `uv` arc. It has a $-1$ entry on the first column entry (indicating leaving the node `u`), and a $+1$ entry on the second column entry (indicating entering the node `v`)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Paths and loops\n", "\n", "Next, we define the terms *paths* and *loops*.\n", "\n", "---\n", "**Definition**\n", "\n", "Given two nodes $x$ and $y$, a *path* from $x$ to $y$ is a sequence $z_{1},z_{2},...,z_{K}$ in $\\mathcal{Z}$ where $z_{1}=x$, $z_{K}=y$, and for every $1\\leq k\\leq K-1$, $\\left( z_{k},z_{k+1}\\right) \\in\\mathcal{A}$.\n", "\n", "\n", "A *loop* (also called `cycle') is a path from a node $x$ to itself.\n", "\n", "---\n", "\n", "**Stylized example**. In our stylized example, `(s,u,v,t)` is a path, but `(s,v,u,t)` is not a path." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Transportation costs\n", "\n", "A vector $c\\in\\mathbb{R}^{\\mathcal{A}}$ defines transportation costs. That is, for $xy\\in\\mathcal{A}$, $c_{xy}$ is the transportation cost associated to arc $xy$. $c$ can also be thought of as the length of arc $xy$. The cost of moving the good from node $x$ to node $y$ along path $z_{1},z_{2},...,z_{K}$ is\n", "\n", "\\begin{align*}\n", "\\sum_{k=1}^{K-1}c_{z_{k}z_{k+1}}.\n", "\\end{align*}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Assumptions about the network\n", "\n", "We shall introduce a number of natural assumptions about the network." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### No profitable loop\n", "---\n", "**Assumption** [No profitable loop]\n", "\n", "\n", "There is no profitable loop, which means that there is no sequence $z_{1},...,z_{K}$ in $\\mathcal{Z}$ such that $z_{K}=z_{1}$, $\\left( z_{k},z_{k+1}\\right) \\in\\mathcal{A}$, and $\\sum_{k=1}^{K-1}c_{z_{k}z_{k+1}}<0$.\n", "\n", "---\n", "\n", "In particular, there is no profitable loop if $c\\geq0$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Balancedness\n", "Let $q_{z}$ be the *net demand* at node $z \\in \\mathcal{Z}$, which is the quantity of goods disappearing from the graph. The set of nodes defined by\n", "\n", "\\begin{align*}\n", "\\mathcal{Z}_{0}=\\left\\{ x\\in\\mathcal{Z}:q_{x}<0\\right\\} \\text{, and }\\mathcal{Z}_{1}=\\left\\{ y\\in\\mathcal{Z}:q_{y}>0\\right\\}\n", "\\end{align*}\n", "\n", "are called the *supply nodes* and *demand nodes* respectively. A node which is neither a supply node, neither a demand node is called a *transit node*.\n", "\n", "Total supply is $-\\sum_{x\\in\\mathcal{Z}_{0}}q_{x}$, total demand is $\\sum_{y\\in\\mathcal{Z}_{1}}q_{y}$.\n", "\n", "---\n", "**Assumption** [Balancedness] \n", "\n", "\n", "Assume that total supply equals total demand on the network, that is\n", "\n", "\\begin{align*}\n", "\\sum_{x\\in\\mathcal{Z}_{0}}q_{x}+\\sum_{y\\in\\mathcal{Z}_{1}}q_{y}=0.\n", "\\end{align*}\n", "\n", "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Connectedness\n", "\n", "---\n", "**Assumption** [Connectedness]\n", "\n", "\n", "Assume the set of supply nodes $\\mathcal{Z}_{0}$ is *strongly connected* to the set of demand nodes $\\mathcal{Z}_{1}$, i.e. for every $x\\in\\mathcal{Z}_{0}$ and $y\\in \\mathcal{Z}_{1}$, there is a path from $x$ to $y$.\n", "\n", "---\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Regular networks\n", "\n", "The specification of the graph, the net demand vector, and the surplus vector\n", "defines a network.\n", "\n", "---\n", "**Definition**\n", "\n", "A directed graph $\\left(\\mathcal{Z},\\mathcal{A}\\right)$ endowed with a net demand vector $\\left(n_{z}\\right)_{z \\in \\mathcal{Z}}$ and a cost vector $\\left(c_{a}\\right)_{a\\in\\mathcal{A}}$ is called a *network* $\\left(\\mathcal{Z}, \\mathcal{A}, q, c\\right)$. If the three assumptions above: [no profitable loop](#noprofitloop), [balancedness](#balancedness) and [connectedness](#connectedness) all hold, the network is called *regular*.\n", "\n", "---\n", "\n", "Without mention of the contrary we shall assume that the network under consideration is regular.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Flows and feasible flows\n", "\n", "A *flow* is a vector $\\mu$ on $\\mathbb{R}_+^{\\mathcal{A}}$, where $\\mu_a$ is the quantity of mass passing through arc $a$.\n", "\n", "The flow of mass disappearing at $z$ equals the flow arriving from other nodes minus the flow shipping to other nodes\n", "\n", "\\begin{align*}\n", "q_{z}=\\sum_{x:xz \\in\\mathcal{A}}\\mu_{xz}-\\sum_{y:zy \\in\\mathcal{A}}\\mu_{zy}\n", "\\end{align*}\n", "\n", "and this equation can be rewritten as $\\nabla^{\\intercal}\\mu=n$. This motivates the following definition:\n", "\n", "---\n", "**Definition**\n", "\n", "The set of feasible flows associated with $q$, denoted $\\mathcal{M}\\left( q\\right)$, or $\\mathcal{M}$ when there is no ambiguity, is defined as the set of flows $\\mu\\geq0$ that verify conservation equation\n", "\n", "\\begin{align*}\n", " \\nabla^{\\intercal}\\mu=q.\n", "\\end{align*}\n", "\n", "---\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Stylized example**. In our stylized example, the flow of mass $\\mu$ associated with arc `(s,u,v,t)` has one unit of mass on arcs `su`,`uv`, and `vt`. Hence, $\\nabla ^\\top \\mu$ is: " ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "array([[-1., 0., 0., 1.]])" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# example_pb.arcsList\n", "# example_pb.nodesList\n", "mu = np.array([1,1,0,0,1])\n", "(example_pb.Nabla().transpose() * mu.reshape(-1,1)).reshape(1,-1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We see that the vector of exiting flows $\\nabla ^\\top \\mu$ is equal to $-1$ on arc `s`, equal to $1$ on arc `t`, and equal to zero elsewhere. Mass ''appears'' at the start of the path `s`, ''disappears'' at the end of the path `t`, and is conserved along the path and elsewhere in the network. \n", "\n", "Conversely, if $q_z = 1\\{z=t\\} - 1\\{z=s\\}$, and for $\\mu\\in\\mathcal{M}(q)$, one may find a path from `s` to `t` along the support of $\\mu$, by the following simple algorithm, which consists of looking for an arc in the support of $\\mu$ originating at the current node, replacing the current node by the destination of the arc, and deleting the arc. It is easy to see that this algorithm will get from `s` to `t` in a finite number of steps. We implement it below: " ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "def pathfinder(self,mu,sind,tind):\n", " # assume mu is such that $\\Nabla^\\top \\mu =1 on t, =-1 on s, and =0 elsewhere\n", " t = self.nodesList[tind]\n", " arcs_active = [(x,y) for (x,y) in self.arcsList if (mu[self.arcsDict[(x,y)]] > 0) ]\n", " current = self.nodesList[sind]\n", " path = [current]\n", " while current != t:\n", " (x,y) = next((x,y) for (x,y) in arcs_active if x==current)\n", " path.append(y)\n", " current = y\n", " arcs_active.remove((x,y))\n", " return(path)\n", " \n", "\n", "network_problem.pathfinder = pathfinder" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Stylized example**. We run this in our stylized example, to determine a path associated with `mu = np.array([1,1,0,0,1])` from `s`to `t`. We get:" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "['s', 'u', 'v', 't']" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "example_pb.pathfinder(mu,0,3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The min-cost flow problem\n", "\n", "Let $p_{z}$ be the price of the commodity at $z$. Consider the strategy which consists in purchasing the good at $x$, shipping to $y$, and selling at $y$. The profit of this strategy is\n", "\n", "\\begin{align*}\n", "p_{y}-p_{x}-c_{xy}=\\left( \\nabla p-c\\right)_{xy}\n", "\\end{align*}\n", "\n", "and hence there is no arbitrage opportunity if $p_{y}-p_{x}-c_{xy}\\leq0$ for every arc $xy$, that is\n", "\n", "\\begin{align*}\n", "\\nabla p\\leq c.\n", "\\end{align*}\n", "\n", "Consider the *minimum cost flow problem*\n", "\n", "\\begin{align*}\n", "& \\min_{\\mu\\geq0}\\sum_{\\left( x,y\\right) \\in\\mathcal{A}}\\mu_{xy}%\n", "c_{xy} \\\\\n", "& s.t.~\\nabla^{\\intercal}\\mu=q\n", "\\end{align*}\n", "\n", "which is a linear programming problem whose dual is\n", "\n", "\\begin{align*}\n", "& \\max_{p\\in\\mathbb{R}^{\\mathcal{Z}}}\\sum_{z\\in\\mathcal{Z}}p\n", "_{z} q_{z}\\\\\n", "& s.t.~\\nabla p\\leq c.\n", "\\end{align*}\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Min-cost flow duality\n", "\n", "---\n", "**Proposition**\n", "1. Under Assumption [No profitable loop](#noprofitloop), the [dual problem](#dual) is feasible, which means that there is a vector $p\\in\\mathbb{R}^{\\mathcal{Z}}$ such that $\\nabla p\\leq c$; and the value of [primal problem](#primal) is strictly less than $+\\infty$.\n", "\n", "2. Under Assumptions [balancedness](#balandedness) and [connectedness](#connectedness), the [primal problem](#primal) is feasible, which means that there is a flow $\\mu\\geq0$ such that $\\nabla^{\\intercal}\\mu=q$; and the value of [dual problem](#dual) is strictly greater than $-\\infty$.\n", "\n", "---\n", "\n", "Assume that $\\left(\\mathcal{Z},\\mathcal{A},q,c\\right)$ is a regular network. Then the value of the [primal problem](#primal) coincides with the value of its [dual](#dual), and both problems have solutions. Further, if $\\mu$ is a solution to the primal and $p$ is a solution to the dual, then $\\mu_{xy}>0$ implies $p_{y}-p_{x}=c_{xy}$." ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "def mincostflow(self,q):\n", " Nabla = self.Nabla()\n", " m=grb.Model()\n", " p = m.addMVar(shape=self.nbNodes, name=\"p\", lb = float('-inf'))\n", " m.setObjective(q @ p, grb.GRB.MAXIMIZE)\n", " m.addConstr( Nabla @ p <= np.array(self.costsList), name=\"Constr\")\n", " m.optimize()\n", " if m.status == grb.GRB.Status.OPTIMAL:\n", " p_z = m.getAttr('x')\n", " mu_a = m.getAttr('pi')\n", " return np.array((p_z,mu_a),dtype=object)\n", "\n", "network_problem.mincostflow = mincostflow" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Special cases\n", "\n", "## Shortest path problem\n", "\n", "Assume there is only one source node $s \\in \\mathcal{Z}$ and one target node $t\\in\\mathcal{Z}$, each associated with unit flow. That is, $q_{z}=1\\left\\{ z=t\\right\\} -1\\left\\{ z=s\\right\\}$. Then the problem boils down to how to push one unit of mass from $s$ to $t$. If we interpret $c_{xy}$ as the distance along arc $xy$, the solution of this problem corresponds to the shortest path from $s$ to $t$. This is why this problem is called shortest path problem. (More generally, this problem extends to the case when $c$ does not have a negative loop).\n", "\n", "The dual problem is then\n", "\n", "\\begin{align*}\n", "\\max_{p} & ~p_{t}-p_{s}\\\\\n", "s.t. & ~p_{y}-p_{x}\\leq c_{xy}~\\forall xy\\in\\mathcal{A}\n", "\\end{align*}\n", "\n", "and we can impose normalization $p_{s}=0$, so that along the travelled path, $p_{z}$ interprets as the distance travelled thus far." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Stylized example**. We run the min-cost flow problem in its ''shortest path'' version on our stylized example, with origin at node `s` and destination at node `t`." ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "scrolled": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Academic license - for non-commercial use only - expires 2022-11-14\n", "Using license file c:\\gurobi911\\gurobi.lic\n", "Gurobi Optimizer version 9.1.1 build v9.1.1rc0 (win64)\n", "Thread count: 6 physical cores, 12 logical processors, using up to 12 threads\n", "Optimize a model with 5 rows, 4 columns and 10 nonzeros\n", "Model fingerprint: 0x930d4ce3\n", "Coefficient statistics:\n", " Matrix range [1e+00, 1e+00]\n", " Objective range [1e+00, 1e+00]\n", " Bounds range [0e+00, 0e+00]\n", " RHS range [1e+00, 4e+00]\n", "Presolve removed 5 rows and 4 columns\n", "Presolve time: 0.00s\n", "Presolve: All rows and columns removed\n", "Iteration Objective Primal Inf. Dual Inf. Time\n", " 0 5.0000000e+00 0.000000e+00 2.000000e-06 0s\n", "\n", "Solved in 0 iterations and 0.01 seconds\n", "Optimal objective 5.000000000e+00\n", "\n", "p_z=[-3.0, 0.0, 0.0, 2.0], and \n", "μ_a =[1.0, 0.0, 0.0, 1.0, -0.0]\n" ] }, { "data": { "text/plain": [ "['s', 'u', 't']" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "p_z,mu_a = example_pb.mincostflow(np.array([[-1., 0., 0., 1.]]))\n", "print('\\np_z='+str(p_z)+', and \\nμ_a ='+str(mu_a) )\n", "#example_pb.arcsList\n", "example_pb.pathfinder(mu_a,0,3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We see that the shortest path in the sytlized example is `sut`, with a total cost of 5." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**NYC subway dataset**. We now move on to the subway dataset, where we are trying to find the shortest path NYU to Columbia." ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Gurobi Optimizer version 9.1.1 build v9.1.1rc0 (win64)\n", "Thread count: 6 physical cores, 12 logical processors, using up to 12 threads\n", "Optimize a model with 1290 rows, 501 columns and 2580 nonzeros\n", "Model fingerprint: 0x15b26335\n", "Coefficient statistics:\n", " Matrix range [1e+00, 1e+00]\n", " Objective range [1e+00, 1e+00]\n", " Bounds range [0e+00, 0e+00]\n", " RHS range [1e-01, 3e+04]\n", "Presolve removed 1131 rows and 309 columns\n", "Presolve time: 0.01s\n", "Presolved: 159 rows, 350 columns, 782 nonzeros\n", "\n", "Iteration Objective Primal Inf. Dual Inf. Time\n", " 0 handle free variables 0s\n", " 148 2.3979662e+04 0.000000e+00 0.000000e+00 0s\n", "\n", "Solved in 148 iterations and 0.01 seconds\n", "Optimal objective 2.397966180e+04\n" ] } ], "source": [ "origin_node = subway_pb.nodesDict['8 St - NYU (N/R)'] \n", "destination_node = subway_pb.nodesDict['116 St - Columbia University (1)'] \n", "q = np.zeros(subway_pb.nbNodes)\n", "q[origin_node] = -1\n", "q[destination_node] = 1\n", "p_z,mu_a = subway_pb.mincostflow(q)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "One can obtain the path by:" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "['8 St - NYU (N/R)',\n", " '14 St - Union Sq (N/R/Q)',\n", " '34 St - Herald Sq (Q/N)',\n", " 'Times Sq - 42 St (R/Q)',\n", " 'Times Sq - 42 St (1/3)',\n", " '72 St (1/2)',\n", " '96 St (2/1/3)',\n", " '103 St (1)',\n", " 'Cathedral Pkwy (1)',\n", " '116 St - Columbia University (1)']" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "subway_pb.pathfinder(mu_a,origin_node,destination_node)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Paris streets dataset**. We shall now look for the shortest path between two emplematic places of the centre of Paris: The Notre Dame cathedral and the Tour Saint-Jacques, through the pedestrian network on the Parisian streets. First, let us use `osmnx`'s function `geocode` to convert the addresses into geographic coordinates, and then `get_nearest_node` to find the two nearest nodes on the network." ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "ename": "NameError", "evalue": "name 'ox' is not defined", "output_type": "error", "traceback": [ "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[1;31mNameError\u001b[0m Traceback (most recent call last)", "\u001b[1;32m\u001b[0m in \u001b[0;36m\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mnyu_geocode\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mox\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgeocode\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'57 boulevard Saint-Germain, 75005 Paris, France'\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;31m# NYU Paris\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 2\u001b[0m \u001b[0mcolumbia_geocode\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mox\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgeocode\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'4 Rue de Chevreuse, 75006 Paris'\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;31m# Columbia in Paris\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 3\u001b[0m \u001b[0morig\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mstreets_pb\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mnodesDict\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mox\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mget_nearest_node\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mstreets_network\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnyu_geocode\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 4\u001b[0m \u001b[0mdest\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mstreets_pb\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mnodesDict\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mox\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mget_nearest_node\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mstreets_network\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcolumbia_geocode\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", "\u001b[1;31mNameError\u001b[0m: name 'ox' is not defined" ] } ], "source": [ "tour_geocode = ox.geocode('Square de la Tour Saint-Jacques, 75004 Paris, France') #Tour Saint-Jacques\n", "notredame_geocode = ox.geocode('6 Parvis Notre-Dame - Pl. Jean-Paul II, 75004 Paris, France') #Notre dame address\n", "orig = streets_pb.nodesDict[ox.get_nearest_node(streets_network, tour_geocode)]\n", "dest = streets_pb.nodesDict[ox.get_nearest_node(streets_network, notredame_geocode)]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We determine the minumum cost flow using Gurobi, and obtain the shortest path accordingly:" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "ename": "NameError", "evalue": "name 'streets_pb' is not defined", "output_type": "error", "traceback": [ "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[1;31mNameError\u001b[0m Traceback (most recent call last)", "\u001b[1;32m\u001b[0m in \u001b[0;36m\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mq\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mzeros\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mstreets_pb\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mnbNodes\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 2\u001b[0m \u001b[0mq\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0morig\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 3\u001b[0m \u001b[0mq\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mdest\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m1\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 4\u001b[0m \u001b[0mp_z\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mmu_a\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mstreets_pb\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmincostflow\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mq\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 5\u001b[0m \u001b[0mpath\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mstreets_pb\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mpathfinder\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mmu_a\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0morig\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mdest\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", "\u001b[1;31mNameError\u001b[0m: name 'streets_pb' is not defined" ] } ], "source": [ "q = np.zeros(streets_pb.nbNodes)\n", "q[orig] = -1\n", "q[dest] = 1\n", "p_z,mu_a = streets_pb.mincostflow(q)\n", "path = streets_pb.pathfinder(mu_a,orig,dest)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can plot the shortest path using `osmnx`'s `plot_graph_route` function." ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "ename": "NameError", "evalue": "name 'ox' is not defined", "output_type": "error", "traceback": [ "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[1;31mNameError\u001b[0m Traceback (most recent call last)", "\u001b[1;32m\u001b[0m in \u001b[0;36m\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mfig\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0max\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mox\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mplot_graph_route\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mstreets_network\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mpath\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mroute_linewidth\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m6\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnode_size\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mbgcolor\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m'k'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m", "\u001b[1;31mNameError\u001b[0m: name 'ox' is not defined" ] } ], "source": [ "fig, ax = ox.plot_graph_route(streets_network, path, route_linewidth=6, node_size=0, bgcolor='k')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Shortest path via dynamic programming\n", "\n", "We have advocated for the use of Gurobi as a black box in this problem, but there exists a direct method to find out the minimal cost path, using dynamic programming.\n", "\n", "The crucial remark here is that if there is a minimal cost path from $s$ to $t$, then there is one of length at most $\\left\\vert \\mathcal{Z}\\right\\vert - 1$.\n", "\n", "For $z\\in\\mathcal{Z}$, and $t\\in\\mathbb{N}$, let $C_{sz}^{k}$ be the minimal cost of the path from $s$ to $z$ among paths of length at most $k$, with the convention that $C_{sz}^{k}=+\\infty$ if there is no such path. One has:\n", "\n", "* $C_{ss}^{0}=0$ and $C_{sz}^{0}=+\\infty$ for all $z\\neq s$, and\n", "\n", "* for $t\\geq1$, $C_{sz}^{k}$ satisfies the Bellman equation:\n", "\n", "\\begin{align*}\n", "C_{sz}^{k}=\\min\\left\\{ C_{sz}^{k},\\min_{x\\in\\mathcal{Z}:sz\\in\\mathcal{A}%\n", "}\\left\\{ c_{sx}+C_{xz}^{k-1}\\right\\} \\right\\} .\n", "\\end{align*}\n", "\n", "It is easy to see that $C_{st}^{\\left\\vert \\mathcal{Z}\\right\\vert }$ is the minimal cost of any path from $s$ to $t$. Shortest-paths algorithms (Dijikstra and Bellman-Ford) implement this idea.\n", "\n", "In `networkx`, this is implemented using:" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "path_bis = nx.shortest_path(streets_network,streets_pb.nodesList[orig], streets_pb.nodesList[dest], weight='length')\n", "path == path_bis" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Exercise**. Translate the above idea into code and recover the shortest path from NYU to Columbia through the NYU subway network." ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [], "source": [ "# your answer here" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Transportation problem\n", "\n", "Assume the problem is bipartite, that is $\\mathcal{Z}=\\mathcal{X}\\cup\\mathcal{Y}$ and $\\mathcal{A}\\subseteq\\mathcal{X}\\times\\mathcal{Y}$. That is, there are no intermediate nodes, and an arc can only go directly from a source to a target.\n", "\n", "Note that any min-cost flow problem can be recast in this form, by taking the shortest distance between any source node and any target node.\n", "\n", "The dual problem is then\n", "\n", "\\begin{align*}\n", "\\max_{p} & ~\\sum_{z\\in\\mathcal{Z}}p_{z}q_{z}\\\\\n", "s.t. & ~p_{y}-p_{x}\\leq c_{xy}~\\forall x\\in\\mathcal{X},y\\in\\mathcal{Y}\n", "\\end{align*}\n", "\n", "which is our first encounter with optimal transport (more on this later in the course)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Soviet planning dataset**. Back to the soviet planning dataset, we can now compute the optimal dispatch of each of the plant." ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Gurobi Optimizer version 9.0.0 build v9.0.0rc2 (win64)\n", "Optimize a model with 155 rows, 78 columns and 310 nonzeros\n", "Model fingerprint: 0x90b94e1d\n", "Coefficient statistics:\n", " Matrix range [1e+00, 1e+00]\n", " Objective range [5e-01, 3e+02]\n", " Bounds range [0e+00, 0e+00]\n", " RHS range [1e+01, 3e+03]\n", "Presolve removed 78 rows and 13 columns\n", "Presolve time: 0.02s\n", "Presolved: 77 rows, 65 columns, 259 nonzeros\n", "\n", "Iteration Objective Primal Inf. Dual Inf. Time\n", " 0 handle free variables 0s\n", " 52 3.9505200e+05 0.000000e+00 0.000000e+00 0s\n", "\n", "Solved in 52 iterations and 0.04 seconds\n", "Optimal objective 3.950520000e+05\n" ] } ], "source": [ "q = np.concatenate((-soviet_df.iloc[68,1:-1],soviet_df.iloc[0:-1,11]))\n", "p_z,mu_a = soviet_pb.mincostflow(q)" ] } ], "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 }