epyt_quality_sensor_placement.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "71a5e074e3fa769b",
   "metadata": {},
   "source": [
    "# Sensor placement for quality monitoring"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2184ea4c6c39bde9",
   "metadata": {},
   "source": [
    "### Initialize EPANET Python Toolkit (EPyT)\n",
    "\n",
    "You should always begin with this command to import the toolkit.\n",
    "\n",
    "[EPyT](https://github.com/OpenWaterAnalytics/EPyT) is available on [PyPI](https://pypi.org/project/epyt/) and can be installed via `pip install epyt`. To upgrade to the latest version if it's already installed, use `pip install --upgrade epyt`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c954362b-a2fb-4127-8e31-f6db33b4026a",
   "metadata": {},
   "outputs": [],
   "source": [
    "%pip install epyt pymoo numpy==1.23"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8bf8122b16425534",
   "metadata": {},
   "outputs": [],
   "source": [
    "from epyt import epanet\n",
    "import numpy as np\n",
    "from math import comb"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c099e86bb8375ee",
   "metadata": {},
   "source": [
    "### Introduction\n",
    "\n",
    "Large-scale networks are difficult to monitor, and it is not possible to have sensors everywhere in the system, mainly due to costs. For this reason, we need to select the best locations to install sensors, in order to guarantee certain criteria, such as:  \n",
    "   - Network coverage  \n",
    "   - Detection delay  \n",
    "   - Impact (e.g., number of people affected)   \n",
    "     \n",
    "                                                      \n",
    "Sensor placement requires the following steps:\n",
    "   1. Create simulation scenarios of contamination events\n",
    "   2. Execute the simulation scenarios\n",
    "   3. Measure the impact\n",
    "   4. Formulate and solve an optimization problem\n",
    "   5. If multiple objectives, costruct the Pareto front and select one solution"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8468433b5305628f",
   "metadata": {},
   "source": [
    "### Load EPANET Network and MSX"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46af525120715f45",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-06-13T09:29:05.622741Z",
     "start_time": "2024-06-13T09:29:05.372396Z"
    }
   },
   "outputs": [],
   "source": [
    "G = epanet('data/BWSN_Network_1.inp')    # Load EPANET Input file"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "500aa64a9a3dad2c",
   "metadata": {},
   "source": [
    "### Simulation Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1538f0f4167c2b24",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-06-13T09:29:05.638302Z",
     "start_time": "2024-06-13T09:29:05.622741Z"
    }
   },
   "outputs": [],
   "source": [
    "t_d = 5 # days\n",
    "G.printv(t_d)\n",
    "\n",
    "G.setTimeSimulationDuration(t_d*24*60*60)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c91eb2ed3e5e967a",
   "metadata": {},
   "source": [
    "### Get Network data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4ecf44dc8ca59dd",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-06-13T09:29:05.665049Z",
     "start_time": "2024-06-13T09:29:05.649819Z"
    }
   },
   "outputs": [],
   "source": [
    "demand_pattern = G.getPattern()\n",
    "roughness_coeff = G.getLinkRoughnessCoeff()\n",
    "node_id = G.getNodeNameID()\n",
    "Nnodes = G.getNodeCount()\n",
    "# *Setup uncertainties*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80af516b70f651f1",
   "metadata": {},
   "source": [
    "### Scenarios"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5564857c9f935e9e",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-06-13T09:29:05.677518Z",
     "start_time": "2024-06-13T09:29:05.665049Z"
    }
   },
   "outputs": [],
   "source": [
    "Ns = 50 # Number of scenarios to simulate\n",
    "u_p = 0.2   # pattern uncertainty\n",
    "u_r = 0.2   # roughness coefficient uncertainty\n",
    "G.printv(u_p)\n",
    "G.printv(u_r)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58b454b5c50cf7d2",
   "metadata": {},
   "source": [
    "##### Create the scenarios - first without contamination, to compute the bounds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d56742e0d79dd9f",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-06-13T09:29:05.910963Z",
     "start_time": "2024-06-13T09:29:05.900213Z"
    }
   },
   "outputs": [],
   "source": [
    "max_inj_conc = 10.0\n",
    "inj_start_time = 1 \n",
    "inj_duration = 24   # maximum duration of 12 hours\n",
    "G.printv(max_inj_conc)\n",
    "G.printv(inj_start_time)\n",
    "G.printv(inj_duration)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f879ba332f748067",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-06-13T09:29:06.392574Z",
     "start_time": "2024-06-13T09:29:06.345165Z"
    }
   },
   "outputs": [],
   "source": [
    "# Generate injection scenarios\n",
    "inj_sc = np.column_stack((\n",
    "    np.random.randint(1, Nnodes + 1, size=Ns),  # Random integers between 1 and Nnodes inclusive\n",
    "    max_inj_conc * np.ones(Ns),\n",
    "    inj_start_time * np.ones(Ns),\n",
    "    inj_duration * np.ones(Ns)\n",
    "))\n",
    "G.printv(inj_sc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a3510bf33e516a8e",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-06-13T09:29:23.476186Z",
     "start_time": "2024-06-13T09:29:06.445524Z"
    }
   },
   "outputs": [],
   "source": [
    "# *Create the coverage matrix*\n",
    "# The coverage matrix is a binary matrix, where earh row is a scenario, and \n",
    "# each column is a node. If a certain contamination scenario affects (can reach) \n",
    "# a certain node, the flag is true, else false.\n",
    "\n",
    "K = np.zeros((Ns, Nnodes))\n",
    "for i in range(Ns):\n",
    "    print(f'Iteration {i+1}')\n",
    "    \n",
    "    # Randomize demands\n",
    "    r_p = -u_p + 2*u_p*np.random.rand(demand_pattern.shape[0], demand_pattern.shape[1])\n",
    "    new_demand_pattern = demand_pattern + demand_pattern*r_p\n",
    "    G.setPatternMatrix(new_demand_pattern) # Set new patterns\n",
    "    \n",
    "    # Randomize pipe roughness\n",
    "    r_r = -u_r + 2*u_r*np.random.rand(len(roughness_coeff))\n",
    "    new_roughness_coeff = roughness_coeff + roughness_coeff*r_r\n",
    "    G.setLinkRoughnessCoeff(new_roughness_coeff) # Set new roughness coefficients\n",
    "    \n",
    "    c = np.zeros(Nnodes)\n",
    "    c[int(inj_sc[i, 0])-1] = inj_sc[i, 1]  \n",
    "    G.setNodeInitialQuality(c.tolist())\n",
    "    \n",
    "    qsp = G.getComputedTimeSeries().NodeQuality\n",
    "    \n",
    "    K[i, :] = np.max(qsp, axis=0) > 0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "266c171a60f00fe4",
   "metadata": {},
   "source": [
    "### Optimization problem\n",
    "  - Avoid having sensors to monitor the same area\n",
    "  - Some scenarios affect all nodes, and some others, may not affect any node.    \n",
    "\n",
    "This means: select those columns that have ones, so that there is at least one positive flag '1' in as many rows as possible  \n",
    "\n",
    "  - $Y   $  is the vector of \"scenario observability\" (number of scenarios × 1)\n",
    "  - $u $  is the solution, a binary vector (number of nodes × 1)\n",
    "  - $K $  is the coverage matrix (number of scenarios × number of nodes) <br>\n",
    "\n",
    "    $ Y = Ku $ \n",
    "\n",
    "The optimization problem is formed (in its simplest form) as  <br>\n",
    "$$ arg max f(u) $$   $$ u \\epsilon \\{0, 1\\}^N $$\n",
    " $$ Σ u_i = m $$    \n",
    "                                                  \n",
    "\n",
    "Where $f(u)$  is the objective function. The simplest function is \"the percentage of scenarios which can be 'detected' by the sensors\".\n",
    "We can have multiple functions, and in that case, we need multi-objective optimization. Solving for a single solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94bcc46b0258b83e",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-06-13T09:29:23.490480Z",
     "start_time": "2024-06-13T09:29:23.481200Z"
    }
   },
   "outputs": [],
   "source": [
    "# random selection of nodes\n",
    "np.random.seed(1)\n",
    "Nsensors = 5    # number of sensors\n",
    "G.printv(Nsensors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2c2468ada13012d",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-06-13T09:29:23.506191Z",
     "start_time": "2024-06-13T09:29:23.493493Z"
    }
   },
   "outputs": [],
   "source": [
    "\n",
    "total_combinations = comb(Nnodes, Nsensors) # 275 Millions for 5 sensors in 129 nodes\n",
    "G.printv(total_combinations)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a309a54901a4926b",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-06-13T09:29:23.539127Z",
     "start_time": "2024-06-13T09:29:23.519180Z"
    }
   },
   "outputs": [],
   "source": [
    "# Create a random solution\n",
    "s = np.random.permutation(Nnodes)[:Nsensors] # return random\n",
    "G.printv(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31ea322f26223532",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-06-13T09:29:23.566117Z",
     "start_time": "2024-06-13T09:29:23.545114Z"
    }
   },
   "outputs": [],
   "source": [
    "# solution\n",
    "u = np.zeros((Nnodes, 1))\n",
    "u[s-1] = 1\n",
    "# compute output\n",
    "Y = np.dot(K, u)\n",
    "G.printv(Y.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6405175cb0701eba",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-06-13T09:29:23.583119Z",
     "start_time": "2024-06-13T09:29:23.570117Z"
    }
   },
   "outputs": [],
   "source": [
    "# How many non-zero elements in Y?\n",
    "perc_coverage = np.sum(Y > 0) / Nsensors * 100\n",
    "G.printv(perc_coverage)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d317ff06aa9b664",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-06-13T09:29:23.598017Z",
     "start_time": "2024-06-13T09:29:23.587118Z"
    }
   },
   "outputs": [],
   "source": [
    "# *Larger scenario dataset*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2219407eb1edfb47",
   "metadata": {},
   "source": [
    "  - How many scenarios should we create?\n",
    "  - In general, the more scenarios, the more robust your solution\n",
    "  - Simulate contamination events at different time steps\n",
    "  - Simulate different hydraulic conditions   \n",
    "  \n",
    "#### Load results from BWSN 2006 competition (37K scenarios)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56eecc96cf0b7daf",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-06-13T09:29:23.679081Z",
     "start_time": "2024-06-13T09:29:23.601005Z"
    }
   },
   "outputs": [],
   "source": [
    "import scipy.io\n",
    "\n",
    "# Load Kall.mat\n",
    "data = scipy.io.loadmat('data/Kall.mat')\n",
    "Kall = data['Kall']\n",
    "\n",
    "# Get the size of Kall\n",
    "Nsall, _ = Kall.shape\n",
    "G.printv(Nsall)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7a2f61849ee4c40",
   "metadata": {},
   "source": [
    "#### Repeat solution computation for the larger dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8dce6234cde2f09c",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-06-13T09:29:23.698580Z",
     "start_time": "2024-06-13T09:29:23.684063Z"
    }
   },
   "outputs": [],
   "source": [
    "Nsensors = 5    # number of sensors\n",
    "G.printv(Nsensors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a67d66f425ae9df8",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-06-13T09:29:23.715115Z",
     "start_time": "2024-06-13T09:29:23.705582Z"
    }
   },
   "outputs": [],
   "source": [
    "# Create a random solution\n",
    "s = np.random.permutation(Nnodes)[:Nsensors]\n",
    "G.printv(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8f3ba3714d5f3a",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-06-13T09:29:23.767367Z",
     "start_time": "2024-06-13T09:29:23.718637Z"
    }
   },
   "outputs": [],
   "source": [
    "# Solution\n",
    "u = np.zeros((Nnodes, 1))\n",
    "u[s - 1] = 1 \n",
    "# Compute output\n",
    "Y = np.dot(Kall, u)\n",
    "# How many non-zero elements in Y?\n",
    "perc_coverage = np.sum(Y > 0) / Nsall * 100\n",
    "G.printv(perc_coverage)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13f8d2539dd4af49",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-06-13T09:29:23.781840Z",
     "start_time": "2024-06-13T09:29:23.772356Z"
    }
   },
   "outputs": [],
   "source": [
    "# *Randomized solutions for multiple sensors*\n",
    "# How many sensors to install? Let's create a lot of random solutions\n",
    " \n",
    "epochs = 1000\n",
    "maxsensors = 20\n",
    "G.printv(epochs)\n",
    "G.printv(maxsensors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "933fbe812515bf1a",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-06-13T09:29:52.181301Z",
     "start_time": "2024-06-13T09:29:23.799880Z"
    }
   },
   "outputs": [],
   "source": [
    "fx = np.zeros((epochs, 2))\n",
    "\n",
    "for i in range(epochs):\n",
    "    print(i)\n",
    "    ns = np.random.randint(1, maxsensors + 1)  \n",
    "    s = np.random.permutation(Nnodes)[:ns]  \n",
    "    u = np.zeros((Nnodes, 1))\n",
    "    u[s - 1] = 1  \n",
    "    Y = np.dot(Kall, u)\n",
    "    fx[i, :] = [ns, np.sum(Y > 0) / Nsall * 100]\n",
    "    \n",
    "print('fx =\\n', fx.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a0f5ca100cba4cac",
   "metadata": {},
   "source": [
    "Plotting the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d4e7423fc1863775",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-06-13T09:29:52.797890Z",
     "start_time": "2024-06-13T09:29:52.181301Z"
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Plot individual points\n",
    "plt.plot(fx[:, 0], fx[:, 1], 'x', markersize=3)\n",
    "\n",
    "# Plot settings\n",
    "plt.grid(True)\n",
    "plt.xlabel('Number of Sensors', fontsize=8)\n",
    "plt.ylabel('Percentage Coverage', fontsize=8)\n",
    "plt.title('Percentage Coverage vs Number of Sensors', fontsize=10)\n",
    "\n",
    "plt.ylim(0, 90)\n",
    "plt.yticks(range(0, 90, 20))\n",
    "\n",
    "plt.xlim(0, 20)\n",
    "plt.xticks(range(0, 21, 4))\n",
    "\n",
    "# Find maximum coverage for each number of sensors\n",
    "maxfx = np.zeros((maxsensors, 2))\n",
    "for i in range(1, maxsensors + 1):\n",
    "    maxfx[i - 1] = [i, np.max(fx[fx[:, 0] == i, 1])]\n",
    "\n",
    "# Plot maximum coverage\n",
    "plt.plot(maxfx[:, 0], maxfx[:, 1], 'r-')\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "daa37de507ae07e9",
   "metadata": {},
   "source": [
    "<b>Question:</b> Is this a valid way of solving problems\n",
    "### Solving the optimization problem\n",
    "  - Integer/Mixed Integer Quadradic Programming\n",
    "  - Evolutionary Optimization (Single-objective)\n",
    "  - Evolutionary Optimization (Multi-objective) Genetic algorithms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6bbdb439715fe1c",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-06-13T09:29:52.812326Z",
     "start_time": "2024-06-13T09:29:52.801521Z"
    }
   },
   "outputs": [],
   "source": [
    "## TESTTTTTTTTTTTTTTTTTTTTTTTTTTT\n",
    "\n",
    "np.random.seed(2)\n",
    "Nsensors = 5\n",
    "print('Nsensors =', Nsensors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db61f497274dbd8a",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-06-13T09:29:52.824954Z",
     "start_time": "2024-06-13T09:29:52.814211Z"
    }
   },
   "outputs": [],
   "source": [
    "Nsall, Nnodes = Kall.shape\n",
    "print('Nsall =', Nsall)\n",
    "print('Nnodes =', Nnodes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12695da21c7f7866",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-06-13T09:29:52.837761Z",
     "start_time": "2024-06-13T09:29:52.825461Z"
    }
   },
   "outputs": [],
   "source": [
    "def gafunc(x):\n",
    "    return calc_score(x, Kall, Nsall)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "953d90e5f83dfa9e",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-06-13T09:29:52.845767Z",
     "start_time": "2024-06-13T09:29:52.838575Z"
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d44d25f788d1b678",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-06-13T09:29:52.856980Z",
     "start_time": "2024-06-13T09:29:52.846091Z"
    }
   },
   "outputs": [],
   "source": [
    "def calc_score(x, Kall, Nsall):\n",
    "    # Calculate the score\n",
    "    scores = np.max(Kall[:, np.round(x).astype(int)], axis=1)\n",
    "    score = -np.sum(scores) / Nsall * 100\n",
    "    return score\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16ce8f69b7dfcacb",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-06-13T09:29:52.870881Z",
     "start_time": "2024-06-13T09:29:52.856980Z"
    }
   },
   "outputs": [],
   "source": [
    "def calc_score_multi(x, Kall, Nsall, Tall):\n",
    "    score1 = -np.sum(np.max(Kall[:, np.round(x).astype(int)], axis=1)) / Nsall * 100\n",
    "    score2 = np.median(np.min(Tall[:, np.round(x).astype(int)], axis=1))\n",
    "    return [score1, score2]"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (Spyder)",
   "language": "python3",
   "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.10.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}