{ "cells": [ { "cell_type": "markdown", "metadata": { "jupyter": { "outputs_hidden": true } }, "source": [ "# Maximal Coverage Location Problem\n", "\n", "*Authors:* [Germano Barcelos](https://github.com/gegen07), [James Gaboardi](https://github.com/jGaboardi), [Levi J. Wolf](https://github.com/ljwolf), [Qunshan Zhao](https://github.com/qszhao)\n", "\n", "LSCP try to minimize the amount of facilities candidate sites in a maximum service standard but then arise another problem: the budget. Sometimes it requires many facilities sites to reach a complete coverage, and there are circumstances when the resources are not available and it's plausible to know how much coverage we can reach using a exact number of facilities. MCLP class try to solve this problem:\n", "\n", "_Maximize the amount of demand covered within a maximal service distance or time standard by locating a fixed number of facilities_\n", "\n", "**MCLP in math notation:**\n", "\n", "$\\begin{array} \\displaystyle \\textbf{Maximize} & \\sum_{i=1}^{n}{a_iy_i} && (1) \\\\\n", "\\displaystyle \\textbf{Subject to:} & \\sum_{j\\in N_i}{x_j \\geq y_i} & \\forall i & (2) \\\\\n", " & \\sum_{j}{x_j = p} & \\forall j & (3) \\\\\n", " & y_i \\in \\{0,1\\} & \\forall i & (4) \\\\\n", " & x_j \\in \\{0,1\\} & \\forall j & (5) \\\\ \\end{array}$\n", "\n", "$\\begin{array} \\displaystyle \\textbf{Where:}\\\\ & & \\displaystyle i & \\small = & \\textrm{index referencing nodes of the network as demand} \\\\\n", "& & j & \\small = & \\textrm{index referencing nodes of the network as potential facility sites} \\\\\n", "& & S & \\small = & \\textrm{maximal acceptable service distance or time standard} \\\\\n", "& & d_{ij} & \\small = & \\textrm{shortest distance or travel time between nodes} i \\textrm{and} j \\\\\n", "& & N_i & \\small = & \\{j | d_{ij} < S\\} \\\\\n", "& & p & \\small = & \\textrm{number of facilities to be located} \\\\\n", "& & x_j & \\small = & \\begin{cases} \n", " 1, \\text{if a facility is located at node } j \\\\\n", " 0, \\text{otherwise} \\\\\n", " \\end{cases} \\\\\n", "& & y_i & \\small = & \\begin{cases} \n", " 1, \\textrm{if demand } i \\textrm{ is covered within a service standard} \\\\\n", " 0, \\textrm{otherwise} \\\\\n", " \\end{cases}\\end{array}$\n", " \n", "_This excerpt above was quoted from Church L., Murray, A. (2018)_\n", "\n", "\n", "This tutorial solves MCLP using spopt.locate.coverage.MCLP instance that depends on a array 2D representing the costs between facilities candidate sites and demand points. For that it uses a lattice 10x10 with simulated points to calculate the costs." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from spopt.locate.coverage import MCLP\n", "from spopt.locate.util import simulated_geo_points\n", "\n", "import numpy\n", "import geopandas\n", "import pulp\n", "import spaghetti\n", "from shapely.geometry import Point\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since the model needs a distance cost matrix we should define some variables. In the comments, it's defined what these variables are for but solver. The solver, assigned below as pulp.PULP_CBC_CMD, is an interface to optimization solver developed by [COIN-OR](https://github.com/coin-or/Cbc). If you want to use another optimization interface as Gurobi or CPLEX see this [guide](https://coin-or.github.io/pulp/guides/how_to_configure_solvers.html) that explains how to achieve this." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "CLIENT_COUNT = 100 # quantity demand points\n", "FACILITY_COUNT = 5 # quantity supply points\n", "\n", "MAX_COVERAGE = 7 # maximum service radius\n", "P_FACILITIES = 4\n", "\n", "# Random seeds for reproducibility\n", "CLIENT_SEED = 5 \n", "FACILITY_SEED = 6 \n", "\n", "solver = pulp.PULP_CBC_CMD(msg=False) # see solvers available in pulp reference" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Lattice 10x10" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create lattice 10x10 with 9 vertical lines in interior." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "lattice = spaghetti.regular_lattice((0, 0, 10, 10), 9, exterior=True)\n", "ntw = spaghetti.Network(in_data=lattice)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Transform spaghetti instance to geopandas geodataframe." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "street = spaghetti.element_as_gdf(ntw, arcs=True)\n", "\n", "street_buffered = geopandas.GeoDataFrame(\n", " geopandas.GeoSeries(street[\"geometry\"].buffer(0.2).unary_union),\n", " crs=street.crs,\n", " columns=[\"geometry\"],\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plotting the network created by spaghetti we can verify that it seems a district with quarters and streets." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD4CAYAAAAq5pAIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAALt0lEQVR4nO3dX4ylBXnH8e+voLFgiZhVQ/nThYTQEhOD3bQoSUNAU9oS6YWNmECIsdleVERrYqFJQ2+aemGMNG0aJ4iSQLDtSgoxRiUoMU3oprNAI7A2GBRYu7pjGqXpjRKfXsyhjusuO3ve9+x5zzzfT0Jmzp95zi+z/OY958w5z6SqkLTz/dKyA0g6NSy71IRll5qw7FITll1q4vRTeWO7du2q3bt3n8qblFo5cODAD6rqDce67JSWfffu3ayvr5/Km5RaSfLc8S7zbrzUhGWXmrDsUhOWXWrCsktNnLDsSe5KciTJk1vOe32Sh5I8M/t49mJjShpqO0f2zwLXHHXercDDVXUx8PDstKQJO+Hv2avq60l2H3X2dcCVs8/vBh4B/nysUO/51KNjjZJ2jH/8k7cN+vp5H7O/qaoOA8w+vvF4V0yyN8l6kvWNjY0TDn7Ppx7l6cMvzhnrFz19+MXR5o05q9u8KWdblXlDD4ILfwVdVa0BawB79uzZ1qaMS885a/BPsZe9/A0aY96Ys7rNm3K2VZo3xLxH9u8nOQdg9vHI4CSSFmresj8I3DT7/CbggXHiSFqU7fzq7T7gUeCSJIeSvB/4GPDOJM8A75ydljRh23k2/r3HuejqkbNIWiBfQSc1YdmlJiy71IRll5qw7FITll1qwrJLTVh2qQnLLjVh2aUmLLvUhGWXmrDsUhOp2tbymFHs2bOnTvS33l5eS3XpOWeNcpsvrwYaY96Ys7rNm3K2VZm3nQ1OSQ5U1Z5jXeaRXWrilP4V1+1yB93OmzflbKs0bwiP7FITll1qwrJLTVh2qQnLLjVh2aUmLLvUhGWXmrDsUhOWXWrCsktNWHapCcsuNWHZpSYsu9TEoLIn+XCSp5I8meS+JK8ZK5ikcc1d9iTnAh8E9lTVm4HTgOvHCiZpXHPvoJuV/d+AtwAvAv8C/G1VfeV4X+MOur7zppxtVeYtbQddVX0X+DjwPHAY+NGxip5kb5L1JOsbGxvz3pykgebeQZfkbOA64ELgh8A/J7mhqu7Zer2qWgPWYPPIvp3Z7qDbefOmnG2V5g0x5Am6dwDfrqqNqvoJcD/w9sGJJC3EkLI/D1ye5IwkAa4GDo4TS9LYhjxm3w/sAx4DvjGbtTZSLkkjG7Q3vqpuB24fKYukBfIVdFITll1qwrJLTVh2qQnLLjVh2aUmLLvUhGWXmrDsUhOWXWrCsktNWHapCcsuNTH3Drp5uIOu77wpZ1uVeUvbQSdptQx6P/uiuINu582bcrZVmjeER3apCcsuNWHZpSYsu9SEZZeasOxSE5ZdasKyS01YdqkJyy41YdmlJiy71IRll5qw7FITg8qe5HVJ9iX5ZpKDScZ5P5+k0Q19P/sdwJeq6t1JXg2cMUImSQsw91qqJGcB/wFcVNsc4lqqvvOmnG1V5i1zLdVFwAbwmSSPJ7kzyZnHuPG9SdaTrG9sbAy4OUlDDLkbfzrwVuDmqtqf5A7gVuAvt16pqtaANdg8sm9nsGupdt68KWdbpXlDDDmyHwIOVdX+2el9bJZf0gTNXfaq+h7wQpJLZmddDTw9SipJoxv6bPzNwL2zZ+KfBd43PJKkRRhU9qp6AjjmM3+SpsVX0ElNWHapCcsuNWHZpSYsu9SEZZeasOxSE5ZdasKyS01YdqkJyy41YdmlJiy71MTcO+jm4Q66vvOmnG1V5i1zB52kFTJ0ecVCuINu582bcrZVmjeER3apCcsuNWHZpSYsu9SEZZeasOxSE5ZdasKyS01YdqkJyy41YdmlJiy71IRll5qw7FITll1qYnDZk5yW5PEkXxgjkKTFGOPIfgtwcIQ5khZo0A66JOcBdwN/DfxZVV37Std3B13feVPOtirzlr2D7pPAR4GfvsKN702ynmR9Y2Nj4M1JmtfcO+iSXAscqaoDSa483vWqag1Yg80j+3Zmu4Nu582bcrZVmjfEkCP7FcC7knwH+BxwVZJ7BieStBBzl72qbquq86pqN3A98NWqumG0ZJJG5e/ZpSZG2RtfVY8Aj4wxS9JieGSXmrDsUhOWXWrCsktNWHapCcsuNWHZpSYsu9SEZZeasOxSE5ZdasKyS01YdqmJQTvoTpY76PrOm3K2VZm37B10klbEKO9nH5s76HbevClnW6V5Q3hkl5qw7FITll1qwrJLTVh2qQnLLjVh2aUmLLvUhGWXmrDsUhOWXWrCsktNWHapCcsuNTF32ZOcn+RrSQ4meSrJLWMGkzSuIe9nfwn4SFU9luRXgANJHqqqp0fKJmlEo62lSvIA8HdV9dDxruNaqr7zppxtVeZNYi1Vkt3AZcD+Y1y2N8l6kvWNjY0xbk7SHAavpUryWuDzwIeq6sWjL6+qNWANNo/s25npWqqdN2/K2VZp3hCDjuxJXsVm0e+tqvsHp5G0MEOejQ/waeBgVX1ivEiSFmHIkf0K4EbgqiRPzP77/ZFySRrZ3I/Zq+pfgYyYRdIC+Qo6qQnLLjVh2aUmLLvUhGWXmrDsUhOWXWrCsktNWHapCcsuNWHZpSYsu9SEZZeaGG0H3Xa4g67vvClnW5V5k9hBJ2n6Bu+gWwR30O28eVPOtkrzhvDILjVh2aUmLLvUhGWXmrDsUhOWXWrCsktNWHapCcsuNWHZpSYsu9SEZZeasOxSE5ZdasKyS00MKnuSa5L8Z5JvJbl1rFCSxjd32ZOcBvw98HvApcB7k1w6VjBJ45p7B12StwF/VVW/Ozt9G0BV/c3xvsYddH3nTTnbqswbuoNuSNnfDVxTVX88O30j8NtV9YGjrrcX2AtwwQUX/OZzzz13wtljrOCRdprtrLh6pbIP2UGXY5z3Cz85qmoNWIPNI/t2Bo+1t0vSzwx5gu4QcP6W0+cB/zUsjqRFGVL2fwcuTnJhklcD1wMPjhNL0tjmvhtfVS8l+QDwZeA04K6qemq0ZJJGNWhvfFV9EfjiSFkkLZCvoJOasOxSE5ZdasKyS02c0j/ZnGQDOPFL6GAX8IMFx5nXlLPBtPNNORvsjHy/VlVvONYFp7Ts25Vk/Xgv+Vu2KWeDaeebcjbY+fm8Gy81YdmlJqZa9rVlB3gFU84G08435Wyww/NN8jG7pPFN9cguaWSWXWpiUmWf8gLLJOcn+VqSg0meSnLLsjMdLclpSR5P8oVlZzlaktcl2Zfkm7Pv4WQ2lCT58Ozf9Mkk9yV5zZLz3JXkSJInt5z3+iQPJXlm9vHsk507mbKvwALLl4CPVNVvAJcDfzqxfAC3AAeXHeI47gC+VFW/DryFieRMci7wQWBPVb2ZzbdrX7/cVHwWuOao824FHq6qi4GHZ6dPymTKDvwW8K2qeraqfgx8DrhuyZn+X1UdrqrHZp//D5v/s5673FQ/k+Q84A+AO5ed5WhJzgJ+B/g0QFX9uKp+uNRQP+904JeTnA6cwZI3LlXV14H/Purs64C7Z5/fDfzhyc6dUtnPBV7YcvoQEyrTVkl2A5cB+5ccZatPAh8FfrrkHMdyEbABfGb2MOPOJGcuOxRAVX0X+DjwPHAY+FFVfWW5qY7pTVV1GDYPPMAbT3bAlMq+rQWWy5bktcDngQ9V1YvLzgOQ5FrgSFUdWHaW4zgdeCvwD1V1GfC/zHE3dBFmj32vAy4EfhU4M8kNy021GFMq++QXWCZ5FZtFv7eq7l92ni2uAN6V5DtsPvy5Ksk9y430cw4Bh6rq5XtC+9gs/xS8A/h2VW1U1U+A+4G3LznTsXw/yTkAs49HTnbAlMo+6QWWScLmY86DVfWJZefZqqpuq6rzqmo3m9+3r1bVZI5OVfU94IUkl8zOuhp4eomRtnoeuDzJGbN/46uZyJOHR3kQuGn2+U3AAyc7YNAOujGtwALLK4AbgW8keWJ23l/M9vDpxG4G7p39IH8WeN+S8wBQVfuT7AMeY/M3Lo+z5JfNJrkPuBLYleQQcDvwMeCfkryfzR9Qf3TSc325rNTDlO7GS1ogyy41YdmlJiy71IRll5qw7FITll1q4v8AXeOi31lykMAAAAAASUVORK5CYII=", "image/svg+xml": [ "\n", "\n", "\n" ], "text/plain": [ "