epyt_contamination_modelling_and_detection.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d0ceeb89-7954-4c5f-8412-c9fd34417bba",
   "metadata": {},
   "source": [
    "## Contamination modelling and detection\n",
    "### Outline\n",
    "\n",
    "- How to simulate a contamination event\n",
    "- Prepare the simulations\n",
    "- Contamination scenario\n",
    "- Simulation the even with MSX\n",
    "- Interval-based contamination detection \n",
    "- Contamination detection"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6172f88a2ab7c778",
   "metadata": {},
   "source": [
    "### Initialize EPANET Python Toolkit (EPyT)\n",
    "#### You should always begin with this command to import the toolkit: `from epyt import epanet`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7402e0b8-c87a-401c-bdbd-2b621930dc5c",
   "metadata": {},
   "outputs": [],
   "source": [
    "%pip install epyt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a3af6b07ca3e90e",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-06-19T08:41:08.892003Z",
     "start_time": "2024-06-19T08:41:02.925001Z"
    }
   },
   "outputs": [],
   "source": [
    "from epyt import epanet\n",
    "import tempfile\n",
    "import requests\n",
    "import os"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8356355ff375694b",
   "metadata": {},
   "source": [
    "### How to simulate a contamination event\n",
    "![Description](icons/img1.png)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "526f4bd114e6aef5",
   "metadata": {},
   "source": [
    "#### Create an arsenite contamination event"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb60a75e-5cca-4eb8-b437-be9ab7b68897",
   "metadata": {},
   "outputs": [],
   "source": [
    "def download_file(url, file_path):\n",
    "    \"\"\"\n",
    "    Downloads a file from the given URL and saves it to the specified file path.\n",
    "\n",
    "    Parameters:\n",
    "    url (str): The URL to download the file from.\n",
    "    file_path (str): The path where the downloaded file should be saved.\n",
    "\n",
    "    Returns:\n",
    "    str: A message indicating the result of the download.\n",
    "    \"\"\"\n",
    "    try:\n",
    "        response = requests.get(url)\n",
    "        response.raise_for_status()  # Check if the request was successful\n",
    "        with open(file_path, 'wb') as file:\n",
    "            file.write(response.content)\n",
    "        return f\"File downloaded and saved as {file_path}\"\n",
    "    except requests.exceptions.HTTPError as err:\n",
    "        return f\"HTTP error occurred: {err}\"\n",
    "    except Exception as err:\n",
    "        return f\"An error occurred: {err}\"\n",
    "\n",
    "# Download BWSN_Network_1.inp and Arsenite.msx\n",
    "url = \"https://raw.githubusercontent.com/WaterFutures/EPyT-and-EPyT-Flow-Tutorial/main/data/BWSN_Network_1.inp\"\n",
    "tmpdirpath = tempfile.gettempdir()\n",
    "inpname = os.path.join(tmpdirpath, \"BWSN_Network_1.inp\")\n",
    "result = download_file(url, inpname)\n",
    "\n",
    "url = \"https://raw.githubusercontent.com/WaterFutures/EPyT-and-EPyT-Flow-Tutorial/main/data/Arsenite.msx\"\n",
    "tmpdirpath = tempfile.gettempdir()\n",
    "msxname = os.path.join(tmpdirpath, \"Arsenite.msx\")\n",
    "resultmsx = download_file(url, msxname)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7deb950f32cab4a",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-06-19T08:41:09.278708Z",
     "start_time": "2024-06-19T08:41:08.892003Z"
    }
   },
   "outputs": [],
   "source": [
    "Gc = epanet(inpname, loadfile=True)   # Load EPANET Input file\n",
    "Gmsx = Gc.loadMSXFile(msxname)        # Load MSX File"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "66a9afd2ca6bcdce",
   "metadata": {},
   "source": [
    "##### Arsenite (As3) reacts with Chlorine (Cl2) and produces Arsenate (As5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ea1a3abc20d7ad1",
   "metadata": {},
   "source": [
    "#### Prepare the simulations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f7715e16a47dcc9",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-06-19T08:41:12.027258Z",
     "start_time": "2024-06-19T08:41:09.288702Z"
    }
   },
   "outputs": [],
   "source": [
    "# Sensor locations\n",
    "sensor_id = ['JUNCTION-17', 'JUNCTION-83', 'JUNCTION-122', 'JUNCTION-31', 'JUNCTION-45']\n",
    "sensor_index = Gc.getNodeIndex(sensor_id)\n",
    "Gc.plot(highlightnode=sensor_index)\n",
    "Gc.plot_close() # Close - duplicate plots in jupyter\n",
    "\n",
    "print(sensor_index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "472f9b61ba12cbe7",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-06-19T08:41:12.042320Z",
     "start_time": "2024-06-19T08:41:12.031373Z"
    }
   },
   "outputs": [],
   "source": [
    "# Simulation Setup\n",
    "t_d = 5 # days\n",
    "Gc.setTimeSimulationDuration(t_d*24*60*60)  #Set simulation duration of 5 days\n",
    "Gc.setMSXTimeStep(3600) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ccf1fd1363abea8",
   "metadata": {},
   "source": [
    "### Contamination scenario\n",
    "#### A contamination scenario needs a location (node), a contaminant concentration and a time profile (start and end time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "931a6a50c2a92b96",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-06-19T08:41:12.059811Z",
     "start_time": "2024-06-19T08:41:12.045363Z"
    }
   },
   "outputs": [],
   "source": [
    "injection_node = 'JUNCTION-23'\n",
    "print('Injection node =', injection_node)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9097be0423de140",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-06-19T08:41:12.073602Z",
     "start_time": "2024-06-19T08:41:12.063819Z"
    }
   },
   "outputs": [],
   "source": [
    "injection_conc = 0.05\n",
    "print('Injection concentration =', injection_conc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee622aa6e1d8fed5",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-06-19T08:41:12.084718Z",
     "start_time": "2024-06-19T08:41:12.076607Z"
    }
   },
   "outputs": [],
   "source": [
    "injection_start_time = 150\n",
    "print('Injection start time =', injection_start_time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dfafa87140d5ade5",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-06-19T08:41:13.096790Z",
     "start_time": "2024-06-19T08:41:12.088723Z"
    }
   },
   "outputs": [],
   "source": [
    "Gc.plot(highlightnode=injection_node);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e5279cc5d40e0399",
   "metadata": {},
   "source": [
    "### Simulation the event with MSX\n",
    "\n",
    "#### Initialize the contamination vector (binary)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d8e2e6dcd201f72",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-06-19T08:41:13.104826Z",
     "start_time": "2024-06-19T08:41:13.096790Z"
    }
   },
   "outputs": [],
   "source": [
    "# Import numpy package\n",
    "import numpy as np\n",
    "\n",
    "# Initialize the vector\n",
    "as3_pat = np.zeros(t_d*48)  \n",
    "\n",
    "# Create the injection pattern \n",
    "as3_pat[injection_start_time:] = 1  \n",
    "\n",
    "print(as3_pat)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8a24b309bc041d5",
   "metadata": {},
   "source": [
    "#### Simulate the contamination of water with arsenite"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bda8b07797fd2b27",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-06-19T09:07:07.567515Z",
     "start_time": "2024-06-19T09:07:07.098429Z"
    }
   },
   "outputs": [],
   "source": [
    "# Specify Arsenite injection source\n",
    "Gc.setMSXSources(injection_node, 'AsIII', 'SETPOINT', injection_conc, 'AS3PAT')\n",
    "\n",
    "# Set pattern for injection\n",
    "Gc.setMSXPattern('AS3PAT',  as3_pat.tolist()) \n",
    "\n",
    "# Solve hydraulics and MSX quality dynamics\n",
    "Qmsx = Gc.getMSXComputedQualityNode(sensor_index)\n",
    "print(Qmsx.Quality[18][:,0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a1a917d972bfd3b2",
   "metadata": {},
   "source": [
    "#### Plot the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9613241ae0d5c837",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "times = Qmsx.Time / 3600\n",
    "for i, index in enumerate(sensor_index):\n",
    "    fig, ax = plt.subplots()\n",
    "    fig.set_figwidth(7)\n",
    "    \n",
    "    plt.plot(times, Qmsx.Quality[index][:,0], label='Cl2 (Chlorine)')\n",
    "    plt.plot(times, Qmsx.Quality[index][:,1], label='As3 (Arsenate)')\n",
    "    plt.plot(times, Qmsx.Quality[index][:,2], label='As5 (Arsenate)')\n",
    "    \n",
    "    plt.legend(loc='upper right', fontsize=7)\n",
    "    plt.title(sensor_id[i], fontsize=8)\n",
    "    plt.xlabel('Time (hours)', fontsize=8)\n",
    "    plt.ylabel('Quality', fontsize=8)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f745aced-87c0-4795-bbae-1955a35b2dc0",
   "metadata": {},
   "source": [
    "##### Unload library"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fbd04e3b-3561-4690-8ae0-005c63e4f512",
   "metadata": {},
   "outputs": [],
   "source": [
    "Gc.unloadMSX()\n",
    "Gc.unload()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ca03eda64220100",
   "metadata": {},
   "source": [
    "### Interval-based contamination detection\n",
    "#### Initialize the simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13039d24caf07038",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load EPANET Network\n",
    "Gc = epanet(inpname) # data/BWSN_Network_1.inp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ab05f39de3b9dd1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load MSX file\n",
    "Gc.loadMSXFile(msxname) # data/Arsenite.msx \n",
    " \n",
    "# Sensor locations\n",
    "sensor_id = ['JUNCTION-17', 'JUNCTION-83', 'JUNCTION-122', 'JUNCTION-31', 'JUNCTION-45']\n",
    "print('sensor_id =', sensor_id)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8d3d5e22597843",
   "metadata": {},
   "outputs": [],
   "source": [
    "sensor_index = Gc.getNodeIndex(sensor_id)\n",
    "print(sensor_index) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d00820f9b7d3654",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set simulation duration for 5 days\n",
    "t_d = 5 # days\n",
    "Gc.setTimeSimulationDuration(t_d*24*60*60)  \n",
    "\n",
    "# Get Network data\n",
    "demand_pattern = Gc.getPattern()\n",
    "roughness_coeff = Gc.getLinkRoughnessCoeff()\n",
    "node_id = Gc.getNodeNameID()\n",
    "\n",
    "Gc.setMSXTimeStep(3600) # in seconds"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b775a375cd32cd2",
   "metadata": {},
   "source": [
    "#### Setup incertainties"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "390d93b3f3890ddf",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-06-19T08:41:14.267827Z",
     "start_time": "2024-06-19T08:41:14.252198Z"
    }
   },
   "outputs": [],
   "source": [
    "# Scenarios\n",
    "Ns = 20 # Number of scenarios to simulate\n",
    "u_p = 0.20  # Pattern uncertainty\n",
    "Gc.printv(u_p) # Print value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "638253fd9ba52bb6",
   "metadata": {},
   "outputs": [],
   "source": [
    "u_r = 0.20  # Roughness coefficient uncertainty\n",
    "Gc.printv(u_r) # Print value"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c028d72f20924f48",
   "metadata": {},
   "source": [
    "#### Create the scenarios - first without contamination, to compute the bounds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96cf80d60f986283",
   "metadata": {},
   "outputs": [],
   "source": [
    "max_inj_conc = 0.0\n",
    "Gc.printv(max_inj_conc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50c54f2abf1a02f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# After day 2 (Dt = 30min)\n",
    "inj_start_time = 2*48   \n",
    "Gc.printv(inj_start_time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae1636be06d94ea1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Maximum duration of 12 hours\n",
    "inj_duration = 24   \n",
    "Gc.printv(inj_duration)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21aee00cd2a6949a",
   "metadata": {},
   "outputs": [],
   "source": [
    "inj_sc = np.column_stack([\n",
    "    np.random.randint(1, Gc.NodeCount + 1, size=Ns),     # Injection locations\n",
    "    max_inj_conc * np.random.rand(Ns),                   # Magnitudes\n",
    "    np.random.randint(1, 49, size=Ns) + inj_start_time,  # Start times\n",
    "    np.random.randint(1, inj_duration + 1, size=Ns)      # Durations\n",
    "])\n",
    "Gc.printv(inj_sc)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a930f92b0c41f367",
   "metadata": {},
   "source": [
    "#### Run epochs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f1d5427ec89b887",
   "metadata": {},
   "outputs": [],
   "source": [
    "QQ = []\n",
    "for i in range(Ns):\n",
    "    print('Iteration', i + 1)\n",
    "    \n",
    "    # Randomize demands\n",
    "    r_p = -u_p + 2 * u_p * np.random.rand(*demand_pattern.shape)\n",
    "    new_demand_pattern = demand_pattern + demand_pattern * r_p\n",
    "    Gc.patterns = new_demand_pattern  # Set new patterns\n",
    "    \n",
    "    # Randomize pipe roughness\n",
    "    r_r = -u_r + 2 * u_r * np.random.rand(*roughness_coeff.shape)\n",
    "    new_roughness_coeff = roughness_coeff + roughness_coeff * r_r\n",
    "    Gc.setLinkRoughnessCoeff(new_roughness_coeff) # Set new roughness coefficients\n",
    "    \n",
    "    # Simulate quality with randomized hydraulic parameters\n",
    "    QQ.append(Gc.getMSXComputedQualityNode(sensor_index))  # Solve hydraulics and MSX quality dynamics\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21e66b8f9e41eac4",
   "metadata": {},
   "source": [
    "#### Plot the results without contaminations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d2819eeeceaa6b8c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create subplots\n",
    "fig, axs = plt.subplots(len(sensor_index), 1, figsize=(10, 15))\n",
    "\n",
    "# Iterate over simulations and sensor indices\n",
    "Qtable = {sensor: [] for sensor in sensor_index}\n",
    "for i in range(Ns):\n",
    "    for j, sensor in enumerate(sensor_index):\n",
    "        # Plot chlorine concentration for each sensor\n",
    "        quality = QQ[i-1].Quality[sensor][:, 0]\n",
    "        axs[j].plot(QQ[i-1].Time / 24 / 60 / 60, quality, '-', color=[0, 0.7, 0.9])\n",
    "        axs[j].grid(True)\n",
    "        axs[j].set_ylabel('Cl_2 (mg/L)')\n",
    "        axs[j].set_xlabel('Time (days)')\n",
    "        Qtable[sensor].append(quality)\n",
    "\n",
    "# Set titles for each subplot\n",
    "for i, ax in enumerate(axs):\n",
    "    ax.set_title(sensor_id[i])\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "print(Qtable)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1aff28cf53311f76",
   "metadata": {},
   "source": [
    "#### Compute the intervals for each sensor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "955a1484-55fd-4c46-b225-0a26a4c7c796",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(Qtable)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "81e825fbef36a347",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize upper and lower bound arrays\n",
    "sensor_len = len(sensor_index)\n",
    "num_time_points = len(Qtable[sensor_index[0]][0])\n",
    "UBQ = np.zeros((num_time_points, sensor_len))\n",
    "LBQ = np.zeros((num_time_points, sensor_len))\n",
    "\n",
    "# Iterate over sensor indices\n",
    "for j, sensor in enumerate(sensor_index):\n",
    "    sensor_data = np.array(Qtable[sensor])\n",
    "\n",
    "    # Calculate upper bound\n",
    "    UBQ[:, j] = np.max(sensor_data, axis=0)\n",
    "    # Calculate lower bound\n",
    "    LBQ[:, j] = np.min(sensor_data, axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "93505769-a7e7-4218-89bf-8fa6b4bcf318",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "# Assuming Qtable has been populated as described earlier\n",
    "\n",
    "# Initialize upper and lower bound arrays\n",
    "sensor_len = len(sensor_index)\n",
    "num_time_points = len(Qtable[sensor_index[0]][0])  # Assuming each quality array has the same length\n",
    "UBQ = np.zeros((num_time_points, sensor_len))\n",
    "LBQ = np.zeros((num_time_points, sensor_len))\n",
    "\n",
    "# Iterate over sensor indices\n",
    "for j, sensor in enumerate(sensor_index):\n",
    "    # Convert list of quality arrays to a numpy array for the current sensor\n",
    "    sensor_data = np.array(Qtable[sensor])\n",
    "    \n",
    "    # Calculate upper bound (max along the time axis)\n",
    "    UBQ[:, j] = np.max(sensor_data, axis=0)\n",
    "    \n",
    "    # Calculate lower bound (min along the time axis)\n",
    "    LBQ[:, j] = np.min(sensor_data, axis=0)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f4eccc8208aae04e",
   "metadata": {},
   "source": [
    "#### Plot the intervals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab93178ab70f340a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create subplots\n",
    "fig, axs = plt.subplots(len(sensor_index), 1, figsize=(10, 15))\n",
    "\n",
    "# Iterate over sensor indices\n",
    "for j, sensor in enumerate(sensor_index):\n",
    "    # Plot upper bound\n",
    "    axs[j].plot(QQ[i].Time / 24 / 60 / 60, UBQ[:, j], 'r-', label='Upper Bound')\n",
    "    # Plot lower bound\n",
    "    axs[j].plot(QQ[i].Time / 24 / 60 / 60, LBQ[:, j], 'r-', label='Lower Bound')\n",
    "    \n",
    "    # Plot the contamination event simulated earlier\n",
    "    axs[j].plot(Qmsx.Time / 24 / 60 / 60, Qmsx.Quality[sensor][:, 0], 'b-', label='Contamination Event')\n",
    "    \n",
    "    axs[j].grid(True)\n",
    "    axs[j].legend(loc='upper right', fontsize=9)\n",
    "    axs[j].set_ylabel('Cl_2 (mg/L)')\n",
    "    axs[j].set_xlabel('Time (days)')\n",
    "    axs[j].set_title(sensor_id[j])\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e34f73068da5405",
   "metadata": {},
   "source": [
    "### Contamination detection\n",
    "##### Compare what the model predicts, with what the (chlorine) sensors measure. If the difference is significantly larger, then put a fault warning.    \n",
    "\n",
    "\n",
    "\n",
    "$ \\\\ \n",
    "e = y - \\hat{y} $\n",
    "\n",
    "If $ e > \\overline{ε} $\n",
    "then fault flag $ φ(k) = 1$\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34c666893e94d6a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create subplots\n",
    "fig, axs = plt.subplots(len(sensor_index), 1, figsize=(10, 15))\n",
    "\n",
    "# Iterate over sensor indices\n",
    "for j, sensor in enumerate(sensor_index):\n",
    "    # Difference between the lower bound and the measured value\n",
    "    e = LBQ[:, j] - Qmsx.Quality[sensor][:, 0]\n",
    "    \n",
    "    # Create the phi signal by putting 1 when the error is > 0\n",
    "    phi = e > 0\n",
    "    \n",
    "    # Plot the results\n",
    "    axs[j].plot(Qmsx.Time / 24 / 60 / 60, phi, 'r-')\n",
    "    axs[j].set_title(sensor_id[j])\n",
    "    axs[j].set_ylabel('φ')\n",
    "    axs[j].set_xlabel('Time (days)')\n",
    "    axs[j].grid(True)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28dce14a-1845-4f38-aaac-8a30896189c3",
   "metadata": {},
   "source": [
    "#### Unload libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c05be48a-e031-496d-9dc5-04b0c792dde6",
   "metadata": {},
   "outputs": [],
   "source": [
    "Gc.unloadMSX()\n",
    "Gc.unload()"
   ]
  }
 ],
 "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": 5
}