epytflow_montecarlo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Uncertainties & Intervals\n",
    "\n",
    "### Monte Carlo Simulations\n",
    "\n",
    "Monte Carlo simulation is a technique used to study how a model  responds to randomly generated inputs. It typically involves a  three-step process: \n",
    "1. Randomly generate “N” inputs (sometimes called scenarios).\n",
    "2. Run a simulation for each of the “N” inputs. Simulations are run on a computerized model of the system being analyzed.\n",
    "3. Aggregate and assess the outputs from the simulations. Common  measures include the mean value of an output, the distribution of output values, and  the minimum or maximum output value.\n",
    "\n",
    "When are Monte Carlo Simulations needed?\n",
    "- Computing *intervals* and *bounds* in signals, when the parameters are uncertain.\n",
    "- *Creating scenarios* for solving various selection problems which cannot be solved by integer or mixed integer programming.\n",
    "\n",
    "### Outline\n",
    "\n",
    "In this notebook we demonstrate how to run (parallel) Monte Carlo Simulations in EPyT-Flow for investigating how uncertainty in the base demands affects the pressures in a network.\n",
    "\n",
    "This is useful for:\n",
    "- Fault diagnosis\n",
    "- State estimation\n",
    "- Forecasting\n",
    "- Sensitivity studies\n",
    "- Demonstrating robustness of solutions\n",
    "- ...\n",
    "\n",
    "We demonstrate how to:\n",
    "1. Setup the scenarios\n",
    "2. Run the Monte Carlo Simulation\n",
    "3. Analyze the results\n",
    "4. Run the Monte Carlo Simulation in parallel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pip install epyt-flow"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\", category=DeprecationWarning)\n",
    "warnings.filterwarnings(\"ignore\", category=ImportWarning)\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from epyt_flow.simulation import ScenarioSimulator, ParallelScenarioSimulation, ScadaData\n",
    "from epyt_flow.uncertainty import ModelUncertainty, UniformUncertainty\n",
    "from epyt_flow.utils import plot_timeseries_data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Setup\n",
    "\n",
    "Prepare the Monte Carlo simulation -- i.e. general parameters and sensor configuration:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Number of simulations\n",
    "n_sim = 100\n",
    "\n",
    "# 5% max uncertainty in base demands\n",
    "eta_bar = 0.05"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We implement the base demand uncertainty by utilizing the EPyT-Flow module on [uncertainties](https://epyt-flow.readthedocs.io/en/stable/tut.uncertainty.html).\n",
    "That is, we derive a new class `MyBaseDemandUncertainty` for implementing the uncertainty logic and (later on) specifying it as part of the [model uncertainty](https://epyt-flow.readthedocs.io/en/stable/epyt_flow.simulation.html#epyt_flow.simulation.scenario_simulator.ScenarioSimulator.set_model_uncertainty):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify and implement the base demand uncertainty\n",
    "# delta = base_demand * uniform_random[-a, a]\n",
    "# base_demand =  base_demand + delta\n",
    "class MyBaseDemandUncertainty(UniformUncertainty):\n",
    "    def __init__(self, **kwds):\n",
    "        super().__init__(**kwds)\n",
    "\n",
    "    def apply(self, data: float) -> float:\n",
    "        z = data * np.random.uniform(low=self.low, high=self.high)\n",
    "        return data + z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify uncertainty\n",
    "base_demand_uncertainty = MyBaseDemandUncertainty(low=-eta_bar, high=eta_bar)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Run the Monte Carlo simulation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run Monte Carlo simulation\n",
    "mcs_results_pressure = []\n",
    "mcs_results_quality = []\n",
    "\n",
    "for i in range(n_sim):\n",
    "    # Create scenario based on Net2\n",
    "    with ScenarioSimulator(f_inp_in=\"Net2.inp\") as sim:\n",
    "        \"\"\"\n",
    "        # TODO: Do it without the model uncertainty class\n",
    "        # Compute and set new base demands\n",
    "        base_demands = sim.epanet_api.getNodeBaseDemands()[1]\n",
    "        delta_bd = (2*np.random.rand(len(base_demands))-1) * eta_bar * base_demands\n",
    "        new_base_demands = base_demands + delta_bd\n",
    "        #print(base_demands)\n",
    "        #print(new_base_demands)\n",
    "\n",
    "        sim.epanet_api.setNodeBaseDemands(new_base_demands)\n",
    "        #\"\"\"\n",
    "        sim.set_model_uncertainty(ModelUncertainty(base_demand_uncertainty=base_demand_uncertainty))\n",
    "\n",
    "        # Place pressure sensors at each node\n",
    "        sim.set_pressure_sensors(sim.sensor_config.nodes)\n",
    "\n",
    "        # Place quality sensors at each node\n",
    "        sim.set_node_quality_sensors(sim.sensor_config.nodes)\n",
    "\n",
    "        # Run simulation and retrieve pressures and quality at each node\n",
    "        scada_data = sim.run_simulation()\n",
    "\n",
    "        #plot_timeseries_data(scada_data.get_data_pressures([\"5\"]).T)\n",
    "        mcs_results_pressure.append(scada_data.get_data_pressures().T)  # Transpose: Each row contains one tim series! \n",
    "        mcs_results_quality.append(scada_data.get_data_nodes_quality().T)\n",
    "\n",
    "# Create NumPy array\n",
    "mcs_results_pressure = np.array(mcs_results_pressure)\n",
    "mcs_results_quality = np.array(mcs_results_quality)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3. Analyze the results\n",
    "\n",
    "How does the pressure fluctuates under the uncertain base demands?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "node_idx = 4   # Investigate the pressure at the fifth node -- refers to node \"5\", recall that indicies start at zero!\n",
    "pressure_at_node = mcs_results_pressure[:, node_idx]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_timeseries_data(pressure_at_node,\n",
    "                     x_axis_label=\"Time steps (1min)\",\n",
    "                     y_axis_label=\"Pressure in $psi$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute upper and lower bounds:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "upper_bound = np.max(pressure_at_node, axis=0)\n",
    "lower_bound = np.min(pressure_at_node, axis=0)\n",
    "average, var = np.mean(pressure_at_node, axis=0), np.var(pressure_at_node, axis=0)\n",
    "\n",
    "upper_bound, lower_bound, average, var"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_, ax = plt.subplots()\n",
    "ax.plot(upper_bound, label=\"Upper bound\")\n",
    "ax.plot(lower_bound, label=\"Lower bound\")\n",
    "ax.plot(average, label=\"Average\")\n",
    "ax.legend()\n",
    "ax.set_xlabel(\"Time steps (1min)\")\n",
    "ax.set_ylabel(\"Pressure in $psi$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Parallel Computations (Advanced topic)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "EPyT-Flow supports the simulation of multiple scenarios in parallel -- see [`ParallelScenarioSimulation`](https://epyt-flow.readthedocs.io/en/stable/epyt_flow.simulation.html#epyt_flow.simulation.parallel_simulation.ParallelScenarioSimulation) for details -- which we are going to utilize to speed up the Monte Carlo simulation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create scenarios:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mcs_sceanrios = []\n",
    "\n",
    "for i in range(n_sim):\n",
    "    # Create sceanrio based on Net2\n",
    "    with ScenarioSimulator(f_inp_in=\"Net2.inp\") as sim:\n",
    "        # Specify base demand uncertainty\n",
    "        base_demand_uncertainty = MyBaseDemandUncertainty(low=-eta_bar, high=eta_bar)\n",
    "        sim.set_model_uncertainty(ModelUncertainty(base_demand_uncertainty=base_demand_uncertainty))\n",
    "\n",
    "        # Specify sensor configuration\n",
    "        sim.set_pressure_sensors(sim.sensor_config.nodes)\n",
    "        sim.set_node_quality_sensors(sim.sensor_config.nodes)\n",
    "\n",
    "        # Export scenario configuration\n",
    "        sim.save_to_epanet_file(f\"Net2_{i}.inp\") # EPyT can not load the same .inp files more than once at the same time!\n",
    "        mcs_sceanrios.append(sim.get_scenario_config())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run Monte Carlo simulation in parallel using up to 4 CPU cores:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Callback handler returns the pressure readings of the node \"5\"\n",
    "def __callback(scada_data: ScadaData, _, scenario_idx: int) -> np.ndarray:\n",
    "    return scada_data.get_data_pressures([\"5\"])\n",
    "\n",
    "# Run simulations in parallel\n",
    "msc_results_pressure = ParallelScenarioSimulation.run(scenarios=mcs_sceanrios, callback=__callback, n_jobs=4)\n",
    "msc_results_pressure = np.array(msc_results_pressure)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Evaluate pressure fluctuations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_timeseries_data(pressure_at_node,\n",
    "                     x_axis_label=\"Time steps (1min)\",\n",
    "                     y_axis_label=\"Pressure in $psi$\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "upper_bound = np.max(pressure_at_node, axis=0)\n",
    "lower_bound = np.min(pressure_at_node, axis=0)\n",
    "average, var = np.mean(pressure_at_node, axis=0), np.var(pressure_at_node, axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "upper_bound, lower_bound, average, var"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "epytflow2",
   "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.9.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}