epyt_joss_advance_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Advanced example\n",
    "\n",
    "A more advanced example is provided below, for designing a simple leakage detection algorithm for the `Net2` benchmark network. The goal is to generate pressure bounds (i.e., the adaptive upper and lower levels of pressure expected at a node, given the uncertainty in model parameters) which can be used to detect events in the system, e.g., by comparing them with available pressure sensor measurements."
   ]
  },
  {
   "cell_type": "markdown",
   "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,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pip install epyt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "from epyt import epanet\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# Create a function to run the simulation and return the pressure results\n",
    "def compute_bounds(G, nsim, base_demands, eta_bar, node_index):\n",
    "    # Seed number to always get the same random results\n",
    "    np.random.seed(1)\n",
    "    # Initialize matrix to save MCS pressures\n",
    "    pmcs = [None for _ in range(nsim)]\n",
    "    for i in range(nsim):\n",
    "        # Compute new base demands\n",
    "        delta_bd = (2 * np.random.rand(1, len(base_demands))[0] - 1) * eta_bar * base_demands\n",
    "        new_base_demands = base_demands + delta_bd\n",
    "        \n",
    "        # Set base demands\n",
    "        G.setNodeBaseDemands(new_base_demands)\n",
    "        \n",
    "        # Compute pressures at each node\n",
    "        pmcs[i] = G.getComputedHydraulicTimeSeries().Pressure\n",
    "        print(f\"Epoch {i}\")\n",
    "\n",
    "    # Compute upper and lower bounds\n",
    "    pmulti = []\n",
    "    for i in range(nsim):\n",
    "        pmulti.append(pmcs[i][:, node_index - 1])\n",
    "    pmulti = np.vstack(pmulti)\n",
    "    ub = np.max(pmulti, axis=0)\n",
    "    lb = np.min(pmulti, axis=0)\n",
    "    meanb = np.mean(pmulti, axis=0)\n",
    "\n",
    "    return pmulti, ub, lb, meanb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For generating leakage events, it’s useful to activate the Pressure-Driven Analysis (PDA), instead of using the default Demand-Driven Analysis (DDA), as the effect on demands due to pressure drops during leakages is not negligible. Moreover, PDA avoids simulation errors due to negative pressures."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "def activate_PDA(G):\n",
    "    type = 'PDA'\n",
    "    pmin = 0\n",
    "    preq = 0.1\n",
    "    pexp = 0.5\n",
    "    G.setDemandModel(type, pmin, preq, pexp)  # Sets the demand model\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "if __name__ == \"__main__\":\n",
    "\n",
    "    # Prepare network for Monte Carlo Simulations\n",
    "    # Load network\n",
    "    inp_name = 'Net2.inp'  # 'L-TOWN.inp'\n",
    "    G = epanet(inp_name)\n",
    "    # Pressure driven analysis\n",
    "    activate_PDA(G)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "    # Get nominal base demands\n",
    "    base_demands = G.getNodeBaseDemands()[1]\n",
    "    print(base_demands)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We assume we have a pressure sensor at the node with ID “11”. We will now create the pressure bounds at that node, using Monte Carlo Simulations (MCS). We assume that there is 2% uncertainty in the nominal base demands compared to the actual demand, which is evenly distributed with the nominal value as the mean. We consider a suitable number of MCS (we use 100 epochs for computational convenience, however, more simulations would provide a more accurate estimation of the bounds). Starting from the current time, we run the simulations for 56 hours for each randomized scenario, the computed pressure measurements are recorded."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "    # Number of simulations\n",
    "    nsim = 100\n",
    "    # Pressure Simulations at Node 5\n",
    "    node_id = '11'\n",
    "    node_index = G.getNodeIndex(node_id)\n",
    "    # 5% max uncertainty in base demands\n",
    "    eta_bar = 0.02\n",
    "    pmulti, ub, lb, meanb = compute_bounds(G, nsim, base_demands, eta_bar, node_index)\n",
    "    print(pmulti, ub, lb, meanb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The upper and lower bounds can be computed by processing all the simulated pressure measurements using numpy methods. The results are depicted in `Figure 3.` Given a sufficient number of simulations, we expect that under normal conditions, pressure at node `“11”` will reside between those bounds. In blue, the average pressure computed by the MCS is depicted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "    # Plots\n",
    "    pressure_units = G.units.NodePressureUnits\n",
    "    plt.rc('xtick', labelsize=7)\n",
    "    plt.rc('ytick', labelsize=7)\n",
    "    fig, ax = plt.subplots(figsize=(4, 3))\n",
    "    ax.plot(ub, 'k')\n",
    "    ax.plot(lb, 'k')\n",
    "    ax.plot(meanb, 'b')\n",
    "    ax.grid(True)\n",
    "    ax.legend(['Upper bound', 'Lower bound', 'Average'], loc='upper right', fontsize=7)\n",
    "    ax.set_title(f'Pressure bounds, Node ID: {node_id}', fontsize=8)\n",
    "    ax.set_xlabel('Time (hours)', fontsize=7)\n",
    "    ax.set_ylabel(f'Pressure ({pressure_units})', fontsize=7)\n",
    "    plt.show()\n",
    "    # fig.savefig('figures/paper_pressure_bounds.png', dpi=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To demonstrate the detection ability of the proposed approach, we simulate a leakage with 50 gallons per minute (GPM) outflow at the node with ID `“7”`, starting `20 hours` after the current time. During a leakage event, we expect that the pressure will drop, and for a sufficiently large leak, the measured pressure can fall below the estimated lower bound, thus triggering a leakage warning.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "    # Add leakage at Node ID 7 after 20 hours\n",
    "    leak_scenario = 50\n",
    "    leak_start = 20\n",
    "    leak_value = 50  # GPM unit\n",
    "    leak_node_id = '7'\n",
    "    leak_node_index = G.getNodeIndex(leak_node_id)\n",
    "    leak_pattern = np.zeros(max(G.getPatternLengths()))\n",
    "    leak_pattern[leak_start:] = 1\n",
    "    pattern_index = G.addPattern('leak', leak_pattern)\n",
    "    G.setNodeDemandPatternIndex(leak_node_index, pattern_index)\n",
    "    G.setNodeBaseDemands(leak_node_index, leak_value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "    # Compute pressures\n",
    "    scada_pressures = G.getComputedHydraulicTimeSeries().Pressure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The detection algorithm compares the lower pressure bound of node `“7”` with the actual pressure as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "    p7 = scada_pressures[:, node_index-1]\n",
    "    e = p7 - lb\n",
    "    alert = e < 0\n",
    "    detectionTime = np.argmax(alert>1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "    # Bounds with Leakage\n",
    "    fig, ax = plt.subplots(figsize=(4, 3))\n",
    "    ax.plot(ub, 'k')\n",
    "    ax.plot(lb, 'k')\n",
    "    ax.plot(p7, 'r')\n",
    "    ax.grid(True)\n",
    "    ax.legend(['Upper bound', 'Lower bound', 'Sensor'], loc='upper right', fontsize=7)\n",
    "    ax.set_title(f'Pressure bounds, Leak Node ID: {leak_node_id}', fontsize=8)\n",
    "    ax.set_xlabel('Time (hours)', fontsize=7)\n",
    "    ax.set_ylabel(f'Pressure ({pressure_units})', fontsize=7)\n",
    "    plt.show()\n",
    "    # fig.savefig('figures/paper_pressure_bounds_leak.png', dpi=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We observe that in this use case, until time `27 hours`, the sensor measurement was within the upper and lower bounds computed in the previous step, therefore there was a `7 hour delay` in detecting the leakage."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "    # Leakage alert\n",
    "    fig, ax = plt.subplots(figsize=(4, 3))\n",
    "    ax.plot(alert)\n",
    "    ax.set_title(f'Leakage alert', fontsize=8)\n",
    "    ax.set_xlabel('Time (hours)', fontsize=7)\n",
    "    plt.show()\n",
    "    # fig.savefig('figures/paper_leakage_alert.png', dpi=300)"
   ]
  }
 ],
 "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"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}