diff --git a/pandapower/estimation/algorithm/base.py b/pandapower/estimation/algorithm/base.py index fd18cce95..cd231680d 100644 --- a/pandapower/estimation/algorithm/base.py +++ b/pandapower/estimation/algorithm/base.py @@ -333,11 +333,15 @@ def estimate(self, eppci: ExtendedPPCI, debug_mode=False, **kwargs): # state vector difference d_E d_E = spsolve(G_m, H.T * (r_inv * r)) + # Scaling of Delta_X to avoid divergence due o ill-conditioning and + # operating conditions far from starting state variables + current_error = np.max(np.abs(d_E)) + if current_error > 0.25: + d_E = d_E*0.25/current_error + # Update E with d_E E += d_E.ravel() - # log data - current_error = np.max(np.abs(d_E)) if debug_mode: obj_func = (r.T*r_inv*r)[0,0] self.logger.debug("Current delta_x: {:.7f}".format(current_error)) diff --git a/tutorials/ukpn_pp_state_estimation.ipynb b/tutorials/ukpn_pp_state_estimation.ipynb new file mode 100644 index 000000000..2f83f69d1 --- /dev/null +++ b/tutorials/ukpn_pp_state_estimation.ipynb @@ -0,0 +1,2708 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "ff5497fc", + "metadata": {}, + "source": [ + "### Pandapower with UK Power Networks - state estimation and forecasting\n", + "\n", + "This tutorial complements the tutorial presented [here](https://github.com/e2nIEE/pandapower/blob/develop/tutorials/ukpn_pp_power_flow.ipynb), which shows how to leverage pandapower for performing studies and analyses on the real grids of UK Power Networks. \n", + "This tutorial shows how to estimate the operating conditions of the grid using ad hoc state estimation techniques able to process measurement or forecast data. \n", + "\n", + "One of the main challenges in distribution systems is to derive the operating conditions of the grid when only a very limited number of measurements is available. \n", + "To this purpose, this tutorial will show the functionalities of the so-called AF-WLS (Allocation Factor based Weighted Least Squares) method available in pandapower. \n", + "This was conceived specifically to deal with scenarios with very few measurements (so-called *unobservable* grids). \n", + "More details about the technical concepts and the mathematical background behind the AF-WLS state estimation algorithm can be found in the publication available at this [link](https://ieeexplore.ieee.org/abstract/document/10497141).\n", + "\n", + "This tutorial has been created in collaboration with UK Power Networks, the Distribution System Operator owning and operating the electricity network across London, the South East and the East of England.\n", + "\n", + "The tutorial will use the real grids associated with the three licensed electricity distribution networks operated by UK Power Networks (LPN, SPN and EPN).\n", + "It will show how to use the AF-WLS state estimation algorithm and its potential performance on reduced portions of the UK Power Networks obtained via the grid reduction algorithm presented in the tutorial available [here](https://github.com/e2nIEE/pandapower/blob/develop/tutorials/ukpn_pp_sensitivity_reduction.ipynb).\n", + "\n", + "UK Power Networks has provided the grid data as part of their LTDS CIM dataset release. It is a \"Shared\" dataset that requires special access. To request access, visit the [LTDS CIM](https://ukpowernetworks.opendatasoft.com/explore/dataset/ukpn-ltds-cim/information/) page and complete the [Shared Data Request Form](https://ukpowernetworks.opendatasoft.com/login/?next=/explore/forms/cim-access-request-form/). Once approved, CIM data is published as XML file attachments (one per licence area: EPN, SPN, LPN). You can download the XML files directly from the portal.\n", + "\n", + "The additional data required to integrate load and generation in the grid are openly available as Excel tables at the following links: \n", + "- EPN --> [EPN Long Term Development Statement - November 2025](https://ukpowernetworks.sharepoint.com/sites/OpenDataPortalLibrary/Shared%20Documents/Forms/AllItems.aspx?id=%2Fsites%2FOpenDataPortalLibrary%2FShared%20Documents%2FGeneral%2FLong%20Term%20Development%20Statement%2FNovember%202025%2FEPN%20Long%20Term%20Development%20Statement%20%2D%20November%202025&p=true&ga=1)\n", + "- SPN --> [SPN Long Term Development Statement - November 2025](https://ukpowernetworks.sharepoint.com/sites/OpenDataPortalLibrary/Shared%20Documents/Forms/AllItems.aspx?id=%2Fsites%2FOpenDataPortalLibrary%2FShared%20Documents%2FGeneral%2FLong%20Term%20Development%20Statement%2FNovember%202025%2FSPN%20Long%20Term%20Development%20Statement%20%2D%20November%202025&p=true&ga=1)\n", + "- LPN --> [LPN Long Term Development Statement - November 2025](https://ukpowernetworks.sharepoint.com/sites/OpenDataPortalLibrary/Shared%20Documents/Forms/AllItems.aspx?id=%2Fsites%2FOpenDataPortalLibrary%2FShared%20Documents%2FGeneral%2FLong%20Term%20Development%20Statement%2FNovember%202025%2FLPN%20Long%20Term%20Development%20Statement%20%2D%20November%202025&p=true&ga=1)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "12e91b69", + "metadata": {}, + "outputs": [], + "source": [ + "# Import the needed libraries \n", + "import pandapower as pp\n", + "import pandapower.topology as top\n", + "from pandapower.toolbox import create_replacement_switch_for_branch, select_subnet\n", + "from sensitivity_functions import build_reduced_network_from_sensitivity\n", + "from pandapower.estimation import state_estimation as se\n", + "from pandapower.estimation.util import remove_shunt_injection_from_meas\n", + "\n", + "import pandas as pd\n", + "import numpy as np\n", + "import copy\n", + "import os\n", + "pd.options.display.float_format = '{:,.4f}'.format\n", + "\n", + "import warnings\n", + "warnings.filterwarnings(\"ignore\")" + ] + }, + { + "cell_type": "markdown", + "id": "cebdfc36", + "metadata": {}, + "source": [ + "#### Workarounds for power flow execution\n", + "The following blocks of code provide some functions to apply some workarounds necessary to run successfully the power flow on the UK Power Networks grids.\n", + "These workarounds include, for example, the creation of external grids (*slack buses* in the power flow terminology) or the replacement of zero impedance components with switches. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "eef635dd", + "metadata": {}, + "outputs": [], + "source": [ + "# Function to replace components with very small impedance with switches.\n", + "from pandapower.toolbox import create_replacement_switch_for_branch\n", + "\n", + "def _replace_zero_impedance_components(net):\n", + " min_ohm = 0.001\n", + " to_replace = (np.abs(net.line.x_ohm_per_km * net.line.length_km) <= min_ohm) & net.line.in_service\n", + "\n", + " if np.any(to_replace):\n", + " print(f\"replaced {sum(to_replace)} lines with switches\")\n", + "\n", + " for i in net.line.loc[to_replace].index.values:\n", + " create_replacement_switch_for_branch(net, \"line\", i)\n", + " net.line.at[i, \"in_service\"] = False\n", + "\n", + " xward = net.xward.loc[(np.abs(net.xward.x_ohm) <= min_ohm) & net.xward.in_service].index.values\n", + " if len(xward) > 0:\n", + " pp.replace_xward_by_ward(net, index=xward, drop=False)\n", + " print(f\"replaced {len(xward)} xwards with wards\")\n", + "\n", + " zb_f_ohm = np.square(net.bus.loc[net.impedance.from_bus.values, \"vn_kv\"].values) / net.impedance.sn_mva\n", + " zb_t_ohm = np.square(net.bus.loc[net.impedance.to_bus.values, \"vn_kv\"].values) / net.impedance.sn_mva\n", + " impedance = ((np.abs(net.impedance.xft_pu) <= min_ohm / zb_f_ohm) |\n", + " (np.abs(net.impedance.xtf_pu) <= min_ohm / zb_t_ohm)) & net.impedance.in_service\n", + "\n", + " if any(impedance):\n", + " print(f\"replaced {sum(impedance)} impedance elements with switches\")\n", + "\n", + " for i in net.impedance.loc[impedance].index.values:\n", + " pp.create_replacement_switch_for_branch(net, \"impedance\", i)\n", + " net.impedance.at[i, \"in_service\"] = False" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "19320e58", + "metadata": {}, + "outputs": [], + "source": [ + "# Function to apply the needed workarounds\n", + "def apply_workarounds(net, license_area, remove_impedance):\n", + " if remove_impedance:\n", + " net.impedance.drop(net.impedance.index, inplace=True)\n", + " _replace_zero_impedance_components(net)\n", + " net.line[\"c_nf_per_km\"] *= 0.1\n", + " net.load[\"p_mw\"] *= 0.1\n", + "\n", + " if license_area == \"LPN\":\n", + " pp.create_ext_grid(net,bus=10711,vm_pu=1)\n", + " pp.create_ext_grid(net,bus=10699,vm_pu=1)\n", + " pp.create_ext_grid(net,bus=10674,vm_pu=1)\n", + " pp.create_ext_grid(net,bus=10738,vm_pu=1)\n", + " pp.create_ext_grid(net,bus=10673,vm_pu=1)\n", + " elif license_area == \"SPN\":\n", + " net.trafo.drop(661,inplace=True)\n", + " pp.create_ext_grid(net,bus=4899,vm_pu=1)\n", + " pp.create_ext_grid(net,bus=4879,vm_pu=1)\n", + " pp.create_ext_grid(net,bus=4903,vm_pu=1)\n", + " pp.create_ext_grid(net,bus=4920,vm_pu=1)\n", + " pp.create_ext_grid(net,bus=4916,vm_pu=1)\n", + " pp.create_ext_grid(net,bus=4925,vm_pu=1)\n", + " pp.create_ext_grid(net,bus=4878,vm_pu=1)\n", + " elif license_area == \"EPN\":\n", + " pp.create_ext_grid(net,bus=9906,vm_pu=1)\n", + " pp.create_ext_grid(net,bus=9918,vm_pu=1)\n", + " pp.create_ext_grid(net,bus=9900,vm_pu=1)\n", + " pp.create_ext_grid(net,bus=9910,vm_pu=1)\n", + " pp.create_ext_grid(net,bus=9878,vm_pu=1)\n", + " else:\n", + " raise ValueError(\"Sorry, this license area does not exist in UK Power Networks. Allowed areas are LPN, SPN and EPN.\")\n", + "\n", + " return net" + ] + }, + { + "cell_type": "markdown", + "id": "38755916", + "metadata": {}, + "source": [ + "#### Sensitivity-based grid reduction\n", + "The following block implements the functions necessary to carry out the grid reduction based on sensitivity factors. \n", + "\n", + "The **goal** of the grid reduction is to reduce the grid around a user-selected bus of interest while keeping, inside the reduced grid, the same power flow behaviour as in the original-size grid. \n", + "\n", + "The main **criterion** for the grid reduction is to cut the grid at transformer level based on the sensitivity of the transformers to the changes applied at the bus of interest. In this way, only the portion of the grid directly affected by changes at the bus of interest is kept within the reduced grid model, whereas other parts of the grid that are not influenced by power variations at the bus of interest are excluded from the model and replaced with equivalent elements. \n", + "\n", + "This grid reduction process allows therefore to create reduced grid models around a selected bus and to focus the analysis on a smaller (and hence more easily manageable) portion of the UK Power Networks grid. \n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "09577b3e", + "metadata": {}, + "outputs": [], + "source": [ + "def calc_trafo_current_sensitivity_from_power_flow(net_start, net_post, min_i_ka=1e-6):\n", + " \"\"\"\n", + " Function to compute the sensitivity of transformers to a power change at the bus of interest\n", + " \"\"\"\n", + " rows = []\n", + " for tidx, tr in net_start.trafo[net_start.trafo.in_service].iterrows():\n", + " tidx = int(tidx)\n", + "\n", + " hv = int(tr.hv_bus)\n", + " lv = int(tr.lv_bus)\n", + "\n", + " # initial currents from PF\n", + " i0_hv_ka_start = abs(float(net_start.res_trafo.i_hv_ka.loc[tidx]))\n", + " i0_lv_ka_start = abs(float(net_start.res_trafo.i_lv_ka.loc[tidx]))\n", + "\n", + " # currents after perturbation from PF\n", + " i0_hv_ka_post = abs(float(net_post.res_trafo.i_hv_ka.loc[tidx]))\n", + " i0_lv_ka_post = abs(float(net_post.res_trafo.i_lv_ka.loc[tidx]))\n", + "\n", + " # current difference between before and after perturbation\n", + " dI_hv_ka = i0_hv_ka_start - i0_hv_ka_post\n", + " dI_lv_ka = i0_lv_ka_start - i0_lv_ka_post\n", + "\n", + " # sensitivity computation\n", + " sf_hv = dI_hv_ka / max(i0_hv_ka_start, float(min_i_ka)) if np.isfinite(dI_hv_ka) else np.nan\n", + " sf_lv = dI_lv_ka / max(i0_lv_ka_start, float(min_i_ka)) if np.isfinite(dI_lv_ka) else np.nan\n", + "\n", + " rows.append({\n", + " \"trafo_index\": tidx,\n", + " \"hv_bus\": hv,\n", + " \"lv_bus\": lv,\n", + " \"vn_hv_kv\": float(net_start.bus.vn_kv.loc[hv]),\n", + " \"vn_lv_kv\": float(net_start.bus.vn_kv.loc[lv]),\n", + " \"i0_hv_ka\": i0_hv_ka_start,\n", + " \"i0_lv_ka\": i0_lv_ka_start,\n", + " \"i0_max_ka\": max(i0_hv_ka_start, i0_lv_ka_start),\n", + " \"dI_hv_ka\": float(dI_hv_ka) if np.isfinite(dI_hv_ka) else np.nan,\n", + " \"dI_lv_ka\": float(dI_lv_ka) if np.isfinite(dI_lv_ka) else np.nan,\n", + " \"dI_max_ka\": max(dI_hv_ka, dI_lv_ka),\n", + " \"sf_hv\": float(sf_hv) if np.isfinite(sf_hv) else np.nan,\n", + " \"sf_lv\": float(sf_lv) if np.isfinite(sf_lv) else np.nan,\n", + " \"sf_max\": float(abs(np.nanmax([sf_hv, sf_lv]))),\n", + " })\n", + "\n", + " return pd.DataFrame(rows).set_index(\"trafo_index\")\n", + "\n", + "\n", + "def calc_trafo3w_current_sensitivity_from_power_flow(net_start, net_post, min_i_ka=1e-6):\n", + " \"\"\"\n", + " Function to compute the sensitivity of 3-winding transformers to a power change at the bus of interest\n", + " \"\"\"\n", + "\n", + " rows = []\n", + " for tidx, tr in net_start.trafo3w[net_start.trafo3w.in_service].iterrows():\n", + " tidx = int(tidx)\n", + "\n", + " hv = int(tr.hv_bus)\n", + " mv = int(tr.mv_bus)\n", + " lv = int(tr.lv_bus)\n", + "\n", + " # initial currents from PF\n", + " i0_hv_ka_start = abs(float(net_start.res_trafo3w.i_hv_ka.loc[tidx])) if \"i_hv_ka\" in net_start.res_trafo3w.columns else np.nan\n", + " i0_mv_ka_start = abs(float(net_start.res_trafo3w.i_mv_ka.loc[tidx])) if \"i_mv_ka\" in net_start.res_trafo3w.columns else np.nan\n", + " i0_lv_ka_start = abs(float(net_start.res_trafo3w.i_lv_ka.loc[tidx])) if \"i_lv_ka\" in net_start.res_trafo3w.columns else np.nan\n", + "\n", + " # currents after perturbation from PF\n", + " i0_hv_ka_post = abs(float(net_post.res_trafo3w.i_hv_ka.loc[tidx])) if \"i_hv_ka\" in net_post.res_trafo3w.columns else np.nan\n", + " i0_mv_ka_post = abs(float(net_post.res_trafo3w.i_mv_ka.loc[tidx])) if \"i_mv_ka\" in net_post.res_trafo3w.columns else np.nan\n", + " i0_lv_ka_post = abs(float(net_post.res_trafo3w.i_lv_ka.loc[tidx])) if \"i_lv_ka\" in net_post.res_trafo3w.columns else np.nan\n", + "\n", + " # current difference between before and after perturbation\n", + " dI_hv_ka = i0_hv_ka_start - i0_hv_ka_post\n", + " dI_mv_ka = i0_mv_ka_start - i0_mv_ka_post\n", + " dI_lv_ka = i0_lv_ka_start - i0_lv_ka_post\n", + "\n", + " # sensitivity computation\n", + " sf_hv = dI_hv_ka / max(i0_hv_ka_start, float(min_i_ka)) if np.isfinite(dI_hv_ka) else np.nan\n", + " sf_mv = dI_mv_ka / max(i0_mv_ka_start, float(min_i_ka)) if np.isfinite(dI_mv_ka) else np.nan\n", + " sf_lv = dI_lv_ka / max(i0_lv_ka_start, float(min_i_ka)) if np.isfinite(dI_lv_ka) else np.nan\n", + "\n", + " rows.append({\n", + " \"trafo3w_index\": tidx,\n", + " \"hv_bus\": hv,\n", + " \"mv_bus\": mv,\n", + " \"lv_bus\": lv,\n", + " \"vn_hv_kv\": float(net_start.bus.vn_kv.loc[hv]),\n", + " \"vn_mv_kv\": float(net_start.bus.vn_kv.loc[mv]),\n", + " \"vn_lv_kv\": float(net_start.bus.vn_kv.loc[lv]),\n", + " \"i0_hv_ka\": float(i0_hv_ka_start) if np.isfinite(i0_hv_ka_start) else np.nan,\n", + " \"i0_mv_ka\": float(i0_mv_ka_start) if np.isfinite(i0_mv_ka_start) else np.nan,\n", + " \"i0_lv_ka\": float(i0_lv_ka_start) if np.isfinite(i0_lv_ka_start) else np.nan,\n", + " \"i0_max_ka\": max(i0_hv_ka_start, i0_mv_ka_start, i0_lv_ka_start),\n", + " \"dI_hv_ka\": float(dI_hv_ka) if np.isfinite(dI_hv_ka) else np.nan,\n", + " \"dI_mv_ka\": float(dI_mv_ka) if np.isfinite(dI_mv_ka) else np.nan,\n", + " \"dI_lv_ka\": float(dI_lv_ka) if np.isfinite(dI_lv_ka) else np.nan,\n", + " \"dI_max_ka\": max(dI_hv_ka, dI_mv_ka, dI_lv_ka),\n", + " \"sf_hv\": float(abs(sf_hv)) if np.isfinite(sf_hv) else np.nan,\n", + " \"sf_mv\": float(abs(sf_mv)) if np.isfinite(sf_mv) else np.nan,\n", + " \"sf_lv\": float(abs(sf_lv)) if np.isfinite(sf_lv) else np.nan,\n", + " \"sf_max\": float(abs(np.nanmax([sf_hv, sf_mv, sf_lv]))),\n", + " })\n", + "\n", + " return pd.DataFrame(rows).set_index(\"trafo3w_index\")\n", + "\n", + "\n", + "def calc_impedance_current_sensitivity_from_power_flow(net_start, net_post, min_i_ka=1e-6):\n", + " \"\"\"\n", + " Function to compute the sensitivity of impedance elements to a power change at the bus of interest. \n", + " Only impedances connecting buses at different voltage levels are taken into account.\n", + " \"\"\"\n", + "\n", + " rows = []\n", + " for iidx, imp in net_start.impedance[net_start.impedance.in_service].iterrows():\n", + " iidx = int(iidx)\n", + " fb = int(imp.from_bus)\n", + " tb = int(imp.to_bus)\n", + "\n", + " fv = net_start.bus.vn_kv.loc[fb]\n", + " tv = net_start.bus.vn_kv.loc[tb]\n", + "\n", + " if fv == tv:\n", + " continue\n", + "\n", + " # initial currents from PF\n", + " i0_from_ka_start = abs(float(net_start.res_impedance.i_from_ka.loc[iidx])) if \"i_from_ka\" in net_start.res_impedance.columns else np.nan\n", + " i0_to_ka_start = abs(float(net_start.res_impedance.i_to_ka.loc[iidx])) if \"i_to_ka\" in net_start.res_impedance.columns else np.nan\n", + "\n", + " # currents after perturbation from PF\n", + " i0_from_ka_post = abs(float(net_post.res_impedance.i_from_ka.loc[iidx])) if \"i_from_ka\" in net_post.res_impedance.columns else np.nan\n", + " i0_to_ka_post = abs(float(net_post.res_impedance.i_to_ka.loc[iidx])) if \"i_to_ka\" in net_post.res_impedance.columns else np.nan\n", + "\n", + " # current difference between before and after perturbation\n", + " dI_from_ka = i0_from_ka_start - i0_from_ka_post\n", + " dI_to_ka = i0_to_ka_start - i0_to_ka_post\n", + "\n", + " # sensitivity computation\n", + " sf_from = dI_from_ka / max(i0_from_ka_start, float(min_i_ka)) if np.isfinite(dI_from_ka) else np.nan\n", + " sf_to = dI_to_ka / max(i0_to_ka_start, float(min_i_ka)) if np.isfinite(dI_to_ka) else np.nan\n", + "\n", + " rows.append({\n", + " \"impedance_index\": iidx,\n", + " \"from_bus\": fb,\n", + " \"to_bus\": tb,\n", + " \"vn_from_kv\": float(net_start.bus.vn_kv.loc[fb]),\n", + " \"vn_to_kv\": float(net_start.bus.vn_kv.loc[tb]),\n", + " \"i0_from_ka\": float(i0_from_ka_start) if np.isfinite(i0_from_ka_start) else np.nan,\n", + " \"i0_to_ka\": float(i0_to_ka_start) if np.isfinite(i0_to_ka_start) else np.nan,\n", + " \"i0_max_ka\": max(i0_from_ka_start, i0_to_ka_start),\n", + " \"dI_from_ka\": float(dI_from_ka) if np.isfinite(dI_from_ka) else np.nan,\n", + " \"dI_to_ka\": float(dI_to_ka) if np.isfinite(dI_to_ka) else np.nan,\n", + " \"dI_max_ka\": max(dI_from_ka, dI_to_ka),\n", + " \"sf_from\": float(abs(sf_from)) if np.isfinite(sf_from) else np.nan,\n", + " \"sf_to\": float(abs(sf_to)) if np.isfinite(sf_to) else np.nan,\n", + " \"sf_max\": float(abs(np.nanmax([sf_from, sf_to]))),\n", + " })\n", + "\n", + " if len(rows):\n", + " return pd.DataFrame(rows).set_index(\"impedance_index\")\n", + "\n", + "\n", + "def create_sets(net):\n", + " \"\"\"\n", + " Function to identify the set of trafo, 3w-trafo and impedance to be considered in the sensitivity analysis.\n", + " \"\"\"\n", + "\n", + " el_pairs = set()\n", + " el_adj = {}\n", + "\n", + " # ------------------------------------------------------------------\n", + " # 2W trafos\n", + " # ------------------------------------------------------------------\n", + " for tidx, tr in net.trafo[net.trafo.in_service].iterrows():\n", + " hv = int(tr.hv_bus)\n", + " lv = int(tr.lv_bus)\n", + "\n", + " el_pairs.add(frozenset((hv, lv)))\n", + " el_adj.setdefault(hv, []).append((lv, (\"trafo\", int(tidx))))\n", + " el_adj.setdefault(lv, []).append((hv, (\"trafo\", int(tidx))))\n", + "\n", + " # ------------------------------------------------------------------\n", + " # 3W trafos\n", + " # ------------------------------------------------------------------\n", + " for tidx, tr in net.trafo3w[net.trafo3w.in_service].iterrows():\n", + " hv = int(tr.hv_bus)\n", + " mv = int(tr.mv_bus)\n", + " lv = int(tr.lv_bus)\n", + "\n", + " # all winding pairs exist electrically\n", + " el_pairs.add(frozenset((hv, mv)))\n", + " el_pairs.add(frozenset((hv, lv)))\n", + " el_pairs.add(frozenset((mv, lv)))\n", + "\n", + " el_adj.setdefault(hv, []).append((mv, (\"trafo3w\", int(tidx))))\n", + " el_adj.setdefault(hv, []).append((lv, (\"trafo3w\", int(tidx))))\n", + " el_adj.setdefault(mv, []).append((hv, (\"trafo3w\", int(tidx))))\n", + " el_adj.setdefault(mv, []).append((lv, (\"trafo3w\", int(tidx))))\n", + " el_adj.setdefault(lv, []).append((hv, (\"trafo3w\", int(tidx))))\n", + " el_adj.setdefault(lv, []).append((mv, (\"trafo3w\", int(tidx))))\n", + "\n", + " # ------------------------------------------------------------------\n", + " # Impedances\n", + " # ------------------------------------------------------------------\n", + " for iidx, imp in net.impedance[net.impedance.in_service].iterrows():\n", + " fb = int(imp.from_bus)\n", + " tb = int(imp.to_bus)\n", + "\n", + " fv = net.bus.vn_kv.loc[fb]\n", + " tv = net.bus.vn_kv.loc[tb]\n", + "\n", + " if fv == tv:\n", + " continue\n", + " \n", + " if fv > tv:\n", + " el_pairs.add(frozenset((fb, tb)))\n", + " else:\n", + " el_pairs.add(frozenset((tb, fb)))\n", + " el_adj.setdefault(fb, []).append((tb, (\"impedance\", int(iidx))))\n", + " el_adj.setdefault(tb, []).append((fb, (\"impedance\", int(iidx))))\n", + "\n", + " return el_pairs, el_adj\n", + "\n", + "\n", + "def cut_by_sensitivity(net, G, start_bus, trafo_sens_df, trafo3w_sens_df=None, impedance_sens_df=None,\n", + " sensitivity_threshold=0.05, min_working_current_ka=1e-4, dI_min_ka=1e-4, vn_max_kv=None, \n", + " respect_switches=True, cut_downward_elements=True, keep_boundary_outside_bus=True):\n", + " \"\"\"\n", + " Traversal from start_bus:\n", + " - elements are cut if:\n", + " sf_max < sensitivity_threshold\n", + " OR dI_max_ka <= dI_min_ka\n", + " OR to_vn > vn_max_kv\n", + "\n", + " Returns\n", + " -------\n", + " kept_buses : set[int]\n", + " boundaries : list[dict]\n", + " \"\"\"\n", + " start_bus = int(start_bus)\n", + " vn = net.bus.vn_kv.astype(float)\n", + "\n", + " el_pairs, el_adj = create_sets(net)\n", + "\n", + " kept_buses = {start_bus}\n", + " visited = {start_bus}\n", + " queue = [start_bus]\n", + " upper_boundaries = []\n", + " lower_boundaries = []\n", + " visited_element_dir = set()\n", + "\n", + " while queue:\n", + " u = queue.pop(0)\n", + "\n", + " # --------------------------------------------------------------\n", + " # 1) non-trafo (or impedance) neighbors\n", + " # --------------------------------------------------------------\n", + " for v in G.neighbors(u):\n", + " v = int(v)\n", + " if (frozenset((u, v)) in el_pairs) or (frozenset((v, u)) in el_pairs):\n", + " continue\n", + " if v not in visited:\n", + " visited.add(v)\n", + " kept_buses.add(v)\n", + " queue.append(v)\n", + "\n", + " # --------------------------------------------------------------\n", + " # 2) trafo or impedance neighbors\n", + " # --------------------------------------------------------------\n", + " for v, el_id in el_adj.get(u, []):\n", + " v = int(v)\n", + " if v > u:\n", + " key = (el_id, int(u), int(v))\n", + " else:\n", + " key = (el_id, int(v), int(u))\n", + " if key in visited_element_dir:\n", + " continue\n", + " visited_element_dir.add(key)\n", + "\n", + " vn_u = float(vn.loc[u])\n", + " vn_v = float(vn.loc[v])\n", + "\n", + " # ==========================================================\n", + " # 2W TRAFO\n", + " # ==========================================================\n", + " if el_id[0] == \"trafo\":\n", + " tidx = int(el_id[1])\n", + "\n", + " # local upward traversal?\n", + " is_upward = vn_v > vn_u + 1e-9\n", + "\n", + " if is_upward or cut_downward_elements:\n", + "\n", + " # sensitivity data\n", + " if tidx in trafo_sens_df.index:\n", + " sf = float(trafo_sens_df.at[tidx, \"sf_max\"]) if \"sf_max\" in trafo_sens_df.columns else np.nan\n", + " i0 = float(trafo_sens_df.at[tidx, \"i0_max_ka\"]) if \"i0_max_ka\" in trafo_sens_df.columns else np.nan\n", + " dI = float(trafo_sens_df.at[tidx, \"dI_max_ka\"]) if \"dI_max_ka\" in trafo_sens_df.columns else np.nan\n", + " else:\n", + " sf = np.nan\n", + " i0 = np.nan\n", + " dI = np.nan\n", + "\n", + " cut_due_to_vn = (\n", + " vn_max_kv is not None \n", + " and vn_v > float(vn_max_kv) + 1e-9\n", + " )\n", + "\n", + " cut_due_to_sens = (\n", + " np.isfinite(sf)\n", + " and np.isfinite(i0)\n", + " and i0 >= float(min_working_current_ka)\n", + " and sf < float(sensitivity_threshold)\n", + " )\n", + "\n", + " cut_due_to_dI_min = (\n", + " np.isfinite(dI)\n", + " and abs(dI) < float(dI_min_ka)\n", + " )\n", + "\n", + " if cut_due_to_vn or cut_due_to_sens or cut_due_to_dI_min:\n", + " tr = net.trafo.loc[tidx]\n", + " trafo_info = {\n", + " \"el_type\": \"trafo\",\n", + " \"el_index\": tidx,\n", + " \"hv_bus\": int(tr.hv_bus),\n", + " \"lv_bus\": int(tr.lv_bus),\n", + " \"boundary_bus_inside\": int(u),\n", + " \"boundary_bus_outside\": int(v),\n", + " \"reason\": \"vn_above_vn_max\" if cut_due_to_vn else \"up_below_current_sensitivity\",\n", + " \"from_vn_kv\": vn_u,\n", + " \"to_vn_kv\": vn_v,\n", + " \"sf_max\": sf,\n", + " \"i0_max_ka\": i0,\n", + " \"dI_max_ka\": dI,\n", + " \"sensitivity_threshold\": float(sensitivity_threshold),\n", + " \"min_working_current_ka\": float(min_working_current_ka),\n", + " \"dI_min_ka\": float(dI_min_ka),\n", + " }\n", + " if is_upward:\n", + " upper_boundaries.append(trafo_info)\n", + " else:\n", + " lower_boundaries.append(trafo_info)\n", + "\n", + " if keep_boundary_outside_bus:\n", + " kept_buses.add(int(v))\n", + " continue\n", + "\n", + " if v not in visited:\n", + " visited.add(v)\n", + " kept_buses.add(v)\n", + " queue.append(v)\n", + "\n", + " continue\n", + "\n", + " # ==========================================================\n", + " # 3W TRAFO\n", + " # ==========================================================\n", + " if el_id[0] == \"trafo3w\":\n", + " tidx = int(el_id[1])\n", + " tr3 = net.trafo3w.loc[tidx]\n", + "\n", + " hv = int(tr3.hv_bus)\n", + " mv = int(tr3.mv_bus)\n", + " lv = int(tr3.lv_bus)\n", + "\n", + " if hv not in {u, v}: \n", + " z = hv\n", + " elif mv not in {u, v}:\n", + " z = mv\n", + " else:\n", + " z = lv\n", + " vn_z = float(vn.loc[z])\n", + "\n", + " if z > u:\n", + " key = (el_id, int(u), int(z))\n", + " else:\n", + " key = (el_id, int(z), int(u))\n", + " visited_element_dir.add(key)\n", + "\n", + " if v > z:\n", + " key = (el_id, int(z), int(v))\n", + " else:\n", + " key = (el_id, int(v), int(z))\n", + " visited_element_dir.add(key)\n", + "\n", + " is_upward = (vn_v > vn_u + 1e-9) or (vn_z > vn_u + 1e-9)\n", + "\n", + " if is_upward or cut_downward_elements:\n", + "\n", + " if trafo3w_sens_df is not None and tidx in trafo3w_sens_df.index:\n", + " sf = float(trafo3w_sens_df.at[tidx, \"sf_max\"]) if \"sf_max\" in trafo3w_sens_df.columns else np.nan\n", + " i0 = float(trafo3w_sens_df.at[tidx, \"i0_max_ka\"]) if \"i0_max_ka\" in trafo3w_sens_df.columns else np.nan\n", + " dI = float(trafo3w_sens_df.at[tidx, \"dI_max_ka\"]) if \"dI_max_ka\" in trafo3w_sens_df.columns else np.nan\n", + " else:\n", + " sf = np.nan\n", + " i0 = np.nan\n", + " dI = np.nan\n", + "\n", + " cut_due_to_vn = (\n", + " vn_max_kv is not None\n", + " and float(vn.loc[hv]) > float(vn_max_kv) + 1e-9\n", + " )\n", + "\n", + " cut_due_to_sens = (\n", + " np.isfinite(sf)\n", + " and np.isfinite(i0)\n", + " and i0 >= float(min_working_current_ka)\n", + " and sf < float(sensitivity_threshold)\n", + " )\n", + "\n", + " cut_due_to_dI_min = (\n", + " np.isfinite(dI)\n", + " and abs(dI) < float(dI_min_ka)\n", + " )\n", + "\n", + " if cut_due_to_vn or cut_due_to_sens or cut_due_to_dI_min:\n", + " trafo_info = {\n", + " \"el_type\": \"trafo3w\",\n", + " \"el_index\": tidx,\n", + " \"hv_bus\": hv,\n", + " \"mv_bus\": mv,\n", + " \"lv_bus\": lv,\n", + " \"boundary_bus_inside\": int(u), \n", + " \"boundary_bus_outside\": int(v),\n", + " \"boundary_bus_other\": int(z),\n", + " \"reason\": \"vn_above_vn_max\" if cut_due_to_vn else \"up_below_current_sensitivity\",\n", + " \"from_vn_kv\": vn_u,\n", + " \"to_vn_kv\": vn_v,\n", + " \"sf_max\": sf,\n", + " \"i0_max_ka\": i0,\n", + " \"dI_max_ka\": dI,\n", + " \"sensitivity_threshold\": float(sensitivity_threshold),\n", + " \"min_working_current_ka\": float(min_working_current_ka),\n", + " \"dI_min_ka\": float(dI_min_ka),\n", + " }\n", + " if is_upward:\n", + " upper_boundaries.append(trafo_info)\n", + " else:\n", + " lower_boundaries.append(trafo_info)\n", + "\n", + " if keep_boundary_outside_bus:\n", + " kept_buses.add(int(v))\n", + " kept_buses.add(int(z))\n", + " continue\n", + "\n", + " if v not in visited:\n", + " visited.add(v)\n", + " visited.add(z)\n", + " kept_buses.add(v)\n", + " kept_buses.add(z)\n", + " queue.append(v)\n", + " queue.append(z)\n", + "\n", + " continue\n", + "\n", + " # ==========================================================\n", + " # IMPEDANCE\n", + " # ==========================================================\n", + " if el_id[0] == \"impedance\":\n", + " iidx = int(el_id[1])\n", + "\n", + " # local upward traversal?\n", + " is_upward = vn_v > vn_u + 1e-9\n", + "\n", + " if is_upward or cut_downward_elements:\n", + "\n", + " # sensitivity data\n", + " if iidx in impedance_sens_df.index:\n", + " sf = float(impedance_sens_df.at[iidx, \"sf_max\"]) if \"sf_max\" in impedance_sens_df.columns else np.nan\n", + " i0 = float(impedance_sens_df.at[iidx, \"i0_max_ka\"]) if \"i0_max_ka\" in impedance_sens_df.columns else np.nan\n", + " dI = float(impedance_sens_df.at[iidx, \"dI_max_ka\"]) if \"dI_max_ka\" in impedance_sens_df.columns else np.nan\n", + " else:\n", + " sf = np.nan\n", + " i0 = np.nan\n", + " dI = np.nan\n", + "\n", + " cut_due_to_vn = (\n", + " vn_max_kv is not None\n", + " and vn_v > float(vn_max_kv) + 1e-9\n", + " )\n", + "\n", + " cut_due_to_sens = (\n", + " np.isfinite(sf)\n", + " and np.isfinite(i0)\n", + " and i0 >= float(min_working_current_ka)\n", + " and sf < float(sensitivity_threshold)\n", + " )\n", + "\n", + " cut_due_to_dI_min = (\n", + " np.isfinite(dI)\n", + " and abs(dI) < float(dI_min_ka)\n", + " )\n", + "\n", + " if cut_due_to_vn or cut_due_to_sens or cut_due_to_dI_min:\n", + " if is_upward:\n", + " hv_bus = v\n", + " lv_bus = u\n", + " else:\n", + " hv_bus = u\n", + " lv_bus = v\n", + "\n", + " imp_info = {\n", + " \"el_type\": \"impedance\",\n", + " \"el_index\": iidx,\n", + " \"hv_bus\": hv_bus,\n", + " \"lv_bus\": lv_bus,\n", + " \"boundary_bus_inside\": int(u),\n", + " \"boundary_bus_outside\": int(v),\n", + " \"reason\": \"vn_above_vn_max\" if cut_due_to_vn else \"up_below_current_sensitivity\",\n", + " \"from_vn_kv\": vn_u,\n", + " \"to_vn_kv\": vn_v,\n", + " \"sf_max\": sf,\n", + " \"i0_max_ka\": i0,\n", + " \"dI_max_ka\": dI,\n", + " \"sensitivity_threshold\": float(sensitivity_threshold),\n", + " \"min_working_current_ka\": float(min_working_current_ka),\n", + " \"dI_min_ka\": float(dI_min_ka),\n", + " }\n", + " if is_upward:\n", + " upper_boundaries.append(imp_info)\n", + " else:\n", + " lower_boundaries.append(imp_info)\n", + "\n", + " continue\n", + "\n", + " if v not in visited:\n", + " visited.add(v)\n", + " kept_buses.add(v)\n", + " queue.append(v)\n", + "\n", + " continue\n", + "\n", + " return kept_buses, upper_boundaries, lower_boundaries\n", + "\n", + "\n", + "def find_trafo_from_ext_grid(net, subnet, G):\n", + " \"\"\"\n", + " This function is used if no external grid exists in the reduced grid.\n", + " It finds the transformer connected to the external grid in the original model.\n", + " \"\"\"\n", + "\n", + " el_pairs, el_adj = create_sets(net)\n", + " \n", + " queue = net.ext_grid[\"bus\"].tolist()\n", + " visited = set(queue)\n", + " visited_element_dir = set()\n", + " created = set()\n", + "\n", + " while queue:\n", + " u = queue.pop(0)\n", + "\n", + " # --------------------------------------------------------------\n", + " # 1) non-trafo neighbors\n", + " # --------------------------------------------------------------\n", + " for v in G.neighbors(u):\n", + " v = int(v)\n", + " if frozenset((u, v)) in el_pairs:\n", + " continue\n", + " if v not in visited:\n", + " visited.add(v)\n", + " queue.append(v)\n", + "\n", + " # --------------------------------------------------------------\n", + " # 2) neighbors\n", + " # --------------------------------------------------------------\n", + " for v, el_id in el_adj.get(u, []):\n", + " v = int(v)\n", + " if v > u:\n", + " key = (el_id, int(u), int(v))\n", + " else:\n", + " key = (el_id, int(v), int(u))\n", + " if key in visited_element_dir:\n", + " continue\n", + " visited_element_dir.add(key)\n", + "\n", + " # ==========================================================\n", + " # 2W TRAFO\n", + " # ==========================================================\n", + " if el_id[0] == \"trafo\":\n", + " tidx = int(el_id[1])\n", + "\n", + " subnet_boundary_trafo = subnet.trafo[subnet.trafo.index == tidx]\n", + " if subnet_boundary_trafo.empty:\n", + " if v not in visited:\n", + " visited.add(v)\n", + " queue.append(v)\n", + " else:\n", + " vm = float(net.res_bus.vm_pu.loc[u])\n", + " va = float(net.res_bus.va_degree.loc[u])\n", + " pp.create_ext_grid(subnet, bus=u, vm_pu=vm, va_degree=va)\n", + " created.add(u)\n", + "\n", + " # ==========================================================\n", + " # 3W TRAFO\n", + " # ==========================================================\n", + " if el_id[0] == \"trafo3w\":\n", + " tidx = int(el_id[1])\n", + "\n", + " subnet_boundary_trafo3w = subnet.trafo3w[subnet.trafo3w.index == tidx]\n", + " if subnet_boundary_trafo3w.empty:\n", + " if v not in visited:\n", + " visited.add(v)\n", + " queue.append(v)\n", + " else:\n", + " vm = float(net.res_bus.vm_pu.loc[u])\n", + " va = float(net.res_bus.va_degree.loc[u])\n", + " pp.create_ext_grid(subnet, bus=u, vm_pu=vm, va_degree=va)\n", + " created.add(u)\n", + "\n", + " # ==========================================================\n", + " # IMPEDANCE\n", + " # ==========================================================\n", + " if el_id[0] == \"impedance\":\n", + " iidx = int(el_id[1])\n", + "\n", + " subnet_boundary_impedance = subnet.impedance[subnet.impedance.index == iidx]\n", + " if subnet_boundary_impedance.empty:\n", + " if v not in visited:\n", + " visited.add(v)\n", + " queue.append(v)\n", + " else:\n", + " vm = float(net.res_bus.vm_pu.loc[u])\n", + " va = float(net.res_bus.va_degree.loc[u])\n", + " pp.create_ext_grid(subnet, bus=u, vm_pu=vm, va_degree=va)\n", + " created.add(u)\n", + "\n", + " return subnet, created\n", + "\n", + "\n", + "def add_boundary_ext_grids(subnet, net, boundaries, subG, G):\n", + " \"\"\"\n", + " Creates ext_grids at the boundary buses that remains inside the subnet.\n", + " It considers the boundary buses met in upstream direction.\n", + " \"\"\"\n", + "\n", + " created = set()\n", + " created_pq = set()\n", + " boundary_list = set()\n", + "\n", + " for bnd in boundaries:\n", + " \n", + " btype = bnd[\"el_type\"]\n", + "\n", + " if (btype == \"trafo3w\") or (btype == \"trafo\"):\n", + " b = int(bnd[\"hv_bus\"])\n", + " else: \n", + " b = int(bnd[\"lv_bus\"])\n", + "\n", + " v = [v for v in subG.neighbors(b)]\n", + " boundary_list.add(b)\n", + "\n", + " vm = float(net.res_bus.vm_pu.loc[b])\n", + " va = float(net.res_bus.va_degree.loc[b])\n", + "\n", + " if b not in created:\n", + " pp.create_ext_grid(subnet, bus=b, vm_pu=vm, va_degree=va)\n", + " created.add(b)\n", + "\n", + " if btype == \"trafo3w\":\n", + " b_pq = bnd[\"boundary_bus_other\"] \n", + " v = [v for v in subG.neighbors(b_pq)] \n", + "\n", + " if len(v) == 2: \n", + " idx = bnd[\"el_index\"]\n", + " if bnd[\"mv_bus\"] == b_pq:\n", + " p = net.res_trafo3w.p_mv_mw.loc[idx]\n", + " q = net.res_trafo3w.q_mv_mvar.loc[idx]\n", + " else:\n", + " p = net.res_trafo3w.p_lv_mw.loc[idx]\n", + " q = net.res_trafo3w.q_lv_mvar.loc[idx]\n", + "\n", + " pp.create_sgen(subnet, bus=b_pq, p_mw=p, q_mvar=q)\n", + " created_pq.add(b_pq)\n", + "\n", + " if subnet.ext_grid.empty:\n", + " subnet, b = find_trafo_from_ext_grid(net, subnet, G)\n", + " created = created.union(b)\n", + "\n", + " return created, created_pq\n", + "\n", + "\n", + "def compensate_pq_inj_for_elements_connected_to_bus(net, pq, b, str):\n", + " \"\"\"\n", + " This function compensates for already existing loads, sgens, or other power injection elements\n", + " already existing at the boundary bus.\n", + " \"\"\"\n", + "\n", + " def compensate_pq(net, b, str, element):\n", + " res_el = \"res_\" + element\n", + "\n", + " if np.any(net[element][net[element].bus==b]):\n", + " if str == \"active\":\n", + " val = net[res_el].p_mw[net[element].bus==b].sum()\n", + " elif str == \"reactive\":\n", + " val = net[res_el].q_mvar[net[element].bus==b].sum()\n", + " else:\n", + " val = 0\n", + "\n", + " return val\n", + "\n", + " pq += compensate_pq(net, b, str, \"load\") # compensation for connected loads\n", + " pq -= compensate_pq(net, b, str, \"sgen\") # compensation for connected sgens\n", + " pq -= compensate_pq(net, b, str, \"gen\") # compensation for connected gens\n", + " pq += compensate_pq(net, b, str, \"shunt\") # compensation for connected shunts\n", + " pq += compensate_pq(net, b, str, \"ward\") # compensation for connected ward\n", + " pq += compensate_pq(net, b, str, \"xward\") # compensation for connected xward\n", + "\n", + " return pq\n", + "\n", + "\n", + "def add_boundary_pq_injections(subnet, net, boundaries, subG):\n", + " \"\"\"\n", + " Creates PQ injections at the boundary bus that remains inside the subnet.\n", + " It considers the boundary buses met in downstream direction.\n", + " \"\"\"\n", + "\n", + " created = set()\n", + "\n", + " for bnd in boundaries:\n", + "\n", + " btype = bnd[\"el_type\"]\n", + " if btype == \"trafo\":\n", + " b = int(bnd[\"lv_bus\"])\n", + " tr_idx = int(bnd[\"el_index\"])\n", + "\n", + " v = [v for v in subG.neighbors(b)]\n", + " create_PQ = True\n", + " if len(v)>1:\n", + " b_volt = subnet.bus.vn_kv.loc[b]\n", + " for it in v:\n", + " v_volt = subnet.bus.vn_kv.loc[it]\n", + " if v_volt == b_volt:\n", + " create_PQ = False\n", + "\n", + " if create_PQ:\n", + " p = float(net.res_trafo.p_lv_mw.loc[tr_idx])\n", + " q = float(net.res_trafo.q_lv_mvar.loc[tr_idx])\n", + "\n", + " if b not in created:\n", + " p = compensate_pq_inj_for_elements_connected_to_bus(net, p, b, \"active\")\n", + " q = compensate_pq_inj_for_elements_connected_to_bus(net, q, b, \"reactive\")\n", + "\n", + " pp.create_sgen(subnet, bus=b, p_mw=p, q_mvar=q)\n", + " created.add(b)\n", + " \n", + " elif bnd[\"el_type\"] == \"trafo3w\":\n", + " b_mv = int(bnd[\"mv_bus\"])\n", + " b_lv = int(bnd[\"lv_bus\"])\n", + " tr_idx = int(bnd[\"el_index\"])\n", + "\n", + " v_mv = [v for v in subG.neighbors(b_mv)]\n", + " if len(v_mv)>2:\n", + " continue\n", + "\n", + " v_lv = [v for v in subG.neighbors(b_lv)]\n", + " if len(v_lv)>2:\n", + " continue\n", + "\n", + " p_mv = float(net.res_trafo3w.p_mv_mw.loc[tr_idx])\n", + " q_mv = float(net.res_trafo3w.q_mv_mvar.loc[tr_idx])\n", + " p_lv = float(net.res_trafo3w.p_lv_mw.loc[tr_idx])\n", + " q_lv = float(net.res_trafo3w.q_lv_mvar.loc[tr_idx])\n", + "\n", + " if b_mv not in created:\n", + " p_mv = compensate_pq_inj_for_elements_connected_to_bus(net, p_mv, b_mv, \"active\")\n", + " q_mv = compensate_pq_inj_for_elements_connected_to_bus(net, q_mv, b_mv, \"reactive\")\n", + "\n", + " if b_lv not in created:\n", + " p_lv = compensate_pq_inj_for_elements_connected_to_bus(net, p_lv, b_lv, \"active\")\n", + " q_lv = compensate_pq_inj_for_elements_connected_to_bus(net, q_lv, b_lv, \"reactive\")\n", + " \n", + " pp.create_sgen(subnet, bus=b_mv, p_mw=p_mv, q_mvar=q_mv)\n", + " pp.create_sgen(subnet, bus=b_lv, p_mw=p_lv, q_mvar=q_lv)\n", + " created.add(b_mv)\n", + " created.add(b_lv)\n", + "\n", + " else:\n", + " b = int(bnd[\"hv_bus\"])\n", + " imp_idx = int(bnd[\"el_index\"])\n", + "\n", + " b_lv = int(bnd[\"lv_bus\"])\n", + " if np.any(subnet.bus.index==b_lv):\n", + " continue\n", + " \n", + " if net.impedance.from_bus.loc[imp_idx] == b:\n", + " p = - float(net.res_impedance.p_from_mw.loc[imp_idx])\n", + " q = - float(net.res_impedance.q_from_mvar.loc[imp_idx])\n", + " else:\n", + " p = - float(net.res_impedance.p_to_mw.loc[imp_idx])\n", + " q = - float(net.res_impedance.q_to_mvar.loc[imp_idx])\n", + "\n", + " pp.create_sgen(subnet, bus=b, p_mw=p, q_mvar=q)\n", + " created.add(b)\n", + "\n", + " return created\n", + "\n", + "\n", + "def build_reduced_network(net, start_bus, method=\"power flow\", sensitivity_threshold=0.05, \n", + " min_working_current_ka=0.01, vn_max_kv=None, deltaP_MW=1.0, deltaQ_Mvar=0.0, \n", + " cut_downward_elements=True):\n", + " \"\"\"\n", + " Complete workflow based on classical bus sensitivity calculation:\n", + " 1) assumes base PF already exists in net\n", + " 2) runs power flows with pwr inj variation\n", + " 3) computes element sensitivities\n", + " 4) cuts elements by sensitivity / vn_max / i_min\n", + " 5) builds subnet\n", + " 6) adds boundary ext_grids or power injections\n", + "\n", + " Returns\n", + " -------\n", + " subnet, trafo_sens_df, trafo3w_sens_df, impedance_sens_df, kept_buses, boundaries, created_ext_grids, created_pq_injections\n", + " \"\"\"\n", + "\n", + " if ~np.any(net.bus.index == start_bus):\n", + " print(\"The selected bus was not found in the considered grid\")\n", + " return net, np.empty, np.empty(0), np.empty(0), np.empty(0), np.empty(0), np.empty(0)\n", + "\n", + " net_post = copy.deepcopy(net)\n", + " pp.runpp(net, run_control=False, max_iteration=100)\n", + "\n", + " pp.create_sgen(net_post, bus=start_bus, p_mw=deltaP_MW, q_mvar=deltaQ_Mvar)\n", + " pp.runpp(net_post, run_control=False, max_iteration=100)\n", + "\n", + " trafo_sens_df = calc_trafo_current_sensitivity_from_power_flow(\n", + " net,\n", + " net_post,\n", + " min_i_ka=1e-6)\n", + " \n", + " trafo3w_sens_df = calc_trafo3w_current_sensitivity_from_power_flow(\n", + " net,\n", + " net_post,\n", + " min_i_ka=1e-6)\n", + " \n", + " impedance_sens_df = calc_impedance_current_sensitivity_from_power_flow(\n", + " net,\n", + " net_post,\n", + " min_i_ka=1e-6)\n", + " \n", + " G = top.create_nxgraph(net, respect_switches=True)\n", + "\n", + " kept_buses, hv_boundaries, lv_boundaries = cut_by_sensitivity(\n", + " net, G,\n", + " start_bus=start_bus,\n", + " trafo_sens_df=trafo_sens_df,\n", + " trafo3w_sens_df=trafo3w_sens_df,\n", + " impedance_sens_df=impedance_sens_df,\n", + " sensitivity_threshold=sensitivity_threshold,\n", + " min_working_current_ka=min_working_current_ka,\n", + " dI_min_ka=1e-4,\n", + " vn_max_kv=vn_max_kv,\n", + " respect_switches=True,\n", + " cut_downward_elements=cut_downward_elements,\n", + " keep_boundary_outside_bus=True)\n", + "\n", + " subnet = select_subnet(net, buses=list(kept_buses), include_results=True)\n", + " subnet.user_pf_options = net.user_pf_options\n", + " subG = top.create_nxgraph(subnet, respect_switches=True)\n", + "\n", + " created_ext_grids, created_pq = add_boundary_ext_grids(subnet, net, hv_boundaries, subG, G)\n", + " created_pq_injections = add_boundary_pq_injections(subnet, net, lv_boundaries, subG)\n", + " created_pq_injections = created_pq_injections.union(created_pq)\n", + "\n", + " boundaries = {}\n", + " boundaries[\"hv\"] = hv_boundaries\n", + " boundaries[\"lv\"] = lv_boundaries\n", + "\n", + " try:\n", + " trafo_sens_sorted = trafo_sens_df.sort_values(by=\"sf_max\", ascending=False)\n", + " except:\n", + " trafo_sens_sorted = None\n", + " try:\n", + " trafo3w_sens_sorted = trafo3w_sens_df.sort_values(by=\"sf_max\", ascending=False)\n", + " except: \n", + " trafo3w_sens_sorted = None\n", + " try:\n", + " impedance_sens_sorted = impedance_sens_df.sort_values(by=\"sf_max\", ascending=False)\n", + " except:\n", + " impedance_sens_sorted = None\n", + "\n", + " return subnet, trafo_sens_sorted, trafo3w_sens_sorted, impedance_sens_sorted, boundaries, created_ext_grids, created_pq_injections" + ] + }, + { + "cell_type": "markdown", + "id": "e952ad0e", + "metadata": {}, + "source": [ + "#### Measurement creation\n", + "The following blocks of code implement the functions necessary to create some measurements in the grid, which are needed for the state estimation or state forecasting purposes. \n", + "\n", + "To emulate real-time state estimation, only a limited number of measurements will be created. \n", + "These will be voltage measurements and active and reactive power measurements at the external grid buses. \n", + "For forecasting purposes, power injection forecast data could be instead available to grid operators. \n", + "\n", + "The following code will implement the funcions to create such measurements as well as the function to add uncertainty to these measurements, which emulates the limited accuracy of measurement instruments (or of forecast data) in the grid. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d2a522c7", + "metadata": {}, + "outputs": [], + "source": [ + "def add_measurement_uncertainty(meas, unc):\n", + " dev = abs(unc/300*meas)\n", + " for m in range(len(meas)):\n", + " meas[m] = np.random.normal(meas[m], dev[m])\n", + " return meas\n", + "\n", + "\n", + "def create_ext_grid_voltage_measurements(net, unc):\n", + " idx = net.ext_grid.bus\n", + " v_values = net.ext_grid[\"vm_pu\"].values\n", + " v_values = add_measurement_uncertainty(v_values, unc)\n", + "\n", + " measv = [\"v\"]*len(v_values)\n", + " element = [\"bus\"]*len(v_values)\n", + " std_dev_v = [unc/300]*len(v_values)\n", + "\n", + " measV = pd.DataFrame({\"measurement_type\":measv, \n", + " \"element_type\":element,\n", + " \"element\":idx,\n", + " \"value\":v_values.tolist(),\n", + " \"std_dev\":std_dev_v})\n", + " \n", + " net[\"measurement\"] = pd.concat([net[\"measurement\"],measV], ignore_index=True)\n", + " return net\n", + "\n", + "\n", + "def create_ext_grid_inj_measurements(net, unc):\n", + " idx = net.ext_grid.bus\n", + " p_values = net.res_bus[\"p_mw\"].loc[net.ext_grid.bus].values\n", + " q_values = net.res_bus[\"q_mvar\"].loc[net.ext_grid.bus].values\n", + " p_values = add_measurement_uncertainty(p_values, unc)\n", + " q_values = add_measurement_uncertainty(q_values, unc)\n", + "\n", + " measp = [\"p\"]*len(p_values)\n", + " measq = [\"q\"]*len(q_values)\n", + " element = [\"bus\"]*len(q_values)\n", + " std_dev_p = [unc/300]*len(p_values)\n", + " std_dev_q = [unc/300]*len(q_values)\n", + "\n", + " measP = pd.DataFrame({\"measurement_type\":measp, \n", + " \"element_type\":element,\n", + " \"element\":idx,\n", + " \"value\":p_values.tolist(),\n", + " \"std_dev\":std_dev_p})\n", + " measQ = pd.DataFrame({\"measurement_type\":measq, \n", + " \"element_type\":element,\n", + " \"element\":idx,\n", + " \"value\":q_values.tolist(),\n", + " \"std_dev\":std_dev_q})\n", + " \n", + " net[\"measurement\"] = pd.concat([net[\"measurement\"],measP,measQ], ignore_index=True)\n", + " return net\n", + "\n", + "\n", + "def create_bus_inj_measurements(net, unc):\n", + " idx = net.bus.index\n", + " p_values = net.res_bus[\"p_mw\"].values\n", + " q_values = net.res_bus[\"q_mvar\"].values\n", + " p_values = add_measurement_uncertainty(p_values, unc)\n", + " q_values = add_measurement_uncertainty(q_values, unc)\n", + "\n", + " measp = [\"p\"]*len(p_values)\n", + " measq = [\"q\"]*len(q_values)\n", + " element = [\"bus\"]*len(q_values)\n", + " std_dev_p = [unc/300]*len(p_values)\n", + " std_dev_q = [unc/300]*len(q_values)\n", + "\n", + " measP = pd.DataFrame({\"measurement_type\":measp, \n", + " \"element_type\":element,\n", + " \"element\":idx,\n", + " \"value\":p_values.tolist(),\n", + " \"std_dev\":std_dev_p})\n", + " measQ = pd.DataFrame({\"measurement_type\":measq, \n", + " \"element_type\":element,\n", + " \"element\":idx,\n", + " \"value\":q_values.tolist(),\n", + " \"std_dev\":std_dev_q})\n", + "\n", + " net[\"measurement\"] = pd.concat([net[\"measurement\"],measP,measQ], ignore_index=True)\n", + " return net" + ] + }, + { + "cell_type": "markdown", + "id": "4b009de7", + "metadata": {}, + "source": [ + "#### UK Power Network grids\n", + "This tutorial assumes that the grids of UK Power Networks have been already imported from the CIM data and saved as pandapower networks in json format. \n", + "To see how to import the UK Power Networks grids starting from the CIM files downloadable from the UK Power Networks portal, please refer to the following [UKPN_CIM2pp_tutorial](). \n", + "Here you can also find how to save the pandapower grid into a json file and how to navigate through the pandapower grid data or the attributes of the different grid components. " + ] + }, + { + "cell_type": "markdown", + "id": "761754a1", + "metadata": {}, + "source": [ + "### State Estimation, Example 1 - LPN grid with very low measurement uncertainty\n", + "\n", + "In this first example, state estimation will be shown on an exemplary portion of the LPN grid. \n", + "\n", + "A very low measurement uncertainty will be applied to test the grid with ideal conditions (measurements very close to true values)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "24e15a1a", + "metadata": {}, + "outputs": [], + "source": [ + "# Import the grid for the analysis\n", + "filename = \"LPN EQ SSH_0401_eq.json\" # Give here the name of the json file with the UKPN grid you want to use\n", + "if os.path.isfile(filename):\n", + " net = pp.from_json(filename)\n", + "else:\n", + " print(\"file does not exist, creating a dummy net\")\n", + " net = pp.create_empty_network()\n", + " bus = pp.create_bus(net, vn_kv=132)\n", + " pp.create_ext_grid(net, bus=bus)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "610fe7fc", + "metadata": {}, + "outputs": [], + "source": [ + "# Apply the workarounds on the selected grid\n", + "license_area = \"LPN\" # Provide here the name of the considered license area. It should be \"LPN\", \"SPN\", or \"EPN\".\n", + "if net.bus.index.size > 1:\n", + " remove_impedance = True # Decide if removing fictious impedances from the grid or not\n", + " net = apply_workarounds(net, license_area, remove_impedance)" + ] + }, + { + "cell_type": "markdown", + "id": "a7b8331b", + "metadata": {}, + "source": [ + "Apply the grid reduction to focus the analysis on a limited portion of the overall network" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "84a43353", + "metadata": {}, + "outputs": [], + "source": [ + "start_bus = 1712 # Select the bus of interest around which you want to reduce the grid\n", + "\n", + "# Call the main function for grid reduction\n", + "subnet, trafo_sens_df, trafo3w_sens_df, impedance_sens_df, boundaries, created_ext_grids, created_pq_injections = build_reduced_network(\n", + " net, \n", + " start_bus=start_bus, # start bus considered for the reduction\n", + " method=\"power_flow\", # used method (only power flow available in this tutorial)\n", + " sensitivity_threshold=0.05, # threshold to decide if cutting or not the subnet\n", + " min_working_current_ka=0.001, # minimum current limit considered for the cutting\n", + " vn_max_kv=50.0, # maximum voltage limit considered for the cutting\n", + " deltaP_MW=1.0, # delta of active power toapplied for the sensitivity calculation\n", + " deltaQ_Mvar=0.0, # delta of reactive power toapplied for the sensitivity calculation\n", + " cut_downward_elements=True) # decide if apply cuts also in downstream direction (lower voltage levels) or not\n", + "\n", + "if ~np.any(net.bus.index == start_bus):\n", + " print(\"The selected bus was not found in the considered grid\")\n", + "else:\n", + " print(\"kept buses:\", len(subnet.bus))\n", + " print(\"boundaries:\", len(boundaries[\"hv\"])+len(boundaries[\"lv\"]))\n", + " print(\"ext_grids created:\", len(created_ext_grids))\n", + " print(\"pq_injections_created:\", len(created_pq_injections))" + ] + }, + { + "cell_type": "markdown", + "id": "c3fa0abf", + "metadata": {}, + "source": [ + "Run a power flow to create the conditions of the grid assumed as reference values and extract the measurements from such conditions." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f86dc7be", + "metadata": {}, + "outputs": [], + "source": [ + "# Run the power flow\n", + "pp.runpp(subnet, run_control=False, lightsim2grid=False, max_iteration=100)\n", + "\n", + "# Clean measurements already present in the grid\n", + "subnet.measurement.drop(subnet.measurement.index, inplace=True)\n", + "\n", + "# Create the voltage measurements with 0.01% uncertainty\n", + "subnet = create_ext_grid_voltage_measurements(subnet, unc=0.01)\n", + "\n", + "# Create the active and reactive power measurements with 0.01% uncertainty\n", + "subnet = create_ext_grid_inj_measurements(subnet, unc=0.01)\n", + "\n", + "# Remove injection values from shunt or ward elements (which would not be seen as power injection measurements)\n", + "remove_shunt_injection_from_meas(subnet,\"shunt\")\n", + "remove_shunt_injection_from_meas(subnet,\"ward\")" + ] + }, + { + "cell_type": "markdown", + "id": "90d44b41", + "metadata": {}, + "source": [ + "Visualize the measurements that will be used for state estimation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3379bb03", + "metadata": {}, + "outputs": [], + "source": [ + "display(subnet.measurement)\n", + "display(\"Total nodes in the grid: \" + str(len(subnet.bus)))\n", + "display(\"Total number of measurements: \" + str(len(subnet.measurement)))\n", + "display(\"Measurement redundancy: \" + \"{:.2f}\".format(100*len(subnet.measurement)/(2*len(subnet.bus))) + \" %\")" + ] + }, + { + "cell_type": "markdown", + "id": "0d61cd13", + "metadata": {}, + "source": [ + "Create load and generation clusters needed for the AF-WLS algorithm.\n", + "\n", + "**Note**: here, for the sake of simplicity, only one cluster for the loads and one cluster for the generators is created. In general, however, the algorithm is able to work with multiple clusters (e.g., residential, commercial, industrial, etc. for loads, or PV, wind, etc. for generation). Those clusters are often available in real grids and are necessary for using the AF-WLS algorithm. In pandapower, those clusters should be assigned within the \"type\" attribute of loads and sgens." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "28f440e3", + "metadata": {}, + "outputs": [], + "source": [ + "# Create the load and generation clusters under the load and sgen \"type\" attribute\n", + "subnet.load[\"type\"] = \"generic_load\"\n", + "subnet.sgen[\"type\"] = \"generic_sgen\"" + ] + }, + { + "cell_type": "markdown", + "id": "50f06fec", + "metadata": {}, + "source": [ + "Modify the nominal values of P and Q for loads and sgens.\n", + "\n", + "**Note**: differently from power flow calculations, state estimation relies only on the use of measurements (provided in the net.measurement table). In the AF-WLS algorithm, the powers of loads and sgens are considered as nominal values, which are adopted to derive the allocation factors estimated in the algorithm. Modifying the values of P and Q, as done in the following block of code, allows therefore to change the \"nominal\" values of power of loads and sgens with respect to those previously considered for creating the reference conditions via the power flow. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "06733034", + "metadata": {}, + "outputs": [], + "source": [ + "# Modify the nominal values of power\n", + "subnet.load[\"p_mw\"] *= 2\n", + "subnet.load[\"q_mvar\"] *= 2\n", + "subnet.sgen[\"p_mw\"] *= 4\n", + "subnet.sgen[\"q_mvar\"] *= 4" + ] + }, + { + "cell_type": "markdown", + "id": "80f505f8", + "metadata": {}, + "source": [ + "Run state estimation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5a1d986a", + "metadata": {}, + "outputs": [], + "source": [ + "# Run the state estimation algorithm\n", + "try: \n", + " success = se.estimate(subnet, algorithm=\"af-wls\", tolerance=1e-4, maximum_iterations=500)\n", + " display(\"State estimation successfully converged in \" + str(success[\"num_iterations\"]) + \" iterations.\")\n", + "except: \n", + " display(\"State estimation did not converge\")" + ] + }, + { + "cell_type": "markdown", + "id": "e546007b", + "metadata": {}, + "source": [ + "Compare voltage magnitude state estimation results to reference results given by the initial power flow:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2942fcb2", + "metadata": {}, + "outputs": [], + "source": [ + "if subnet.res_bus_est.index.size > 1: \n", + " volt_diff = subnet.res_bus[\"vm_pu\"].values - subnet.res_bus_est[\"vm_pu\"].values\n", + " max_volt_diff = np.max(abs(volt_diff))\n", + " display(\"Maximum voltage difference = \" + \"{:.6f}\".format(max_volt_diff) + \" p.u.\")" + ] + }, + { + "cell_type": "markdown", + "id": "f005198b", + "metadata": {}, + "source": [ + "Compare line current magnitude state estimation results to reference results given by the initial power flow:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1cff07e5", + "metadata": {}, + "outputs": [], + "source": [ + "if subnet.res_line_est.index.size > 1: \n", + " curr_diff = subnet.res_line[\"i_from_ka\"].values - subnet.res_line_est[\"i_from_ka\"].values\n", + " max_curr_diff = np.max(abs(curr_diff))\n", + " display(\"Maximum current difference = \" + \"{:.6f}\".format(max_curr_diff) + \" kA\")" + ] + }, + { + "cell_type": "markdown", + "id": "3e9036fa", + "metadata": {}, + "source": [ + "Visualize estimated allocation factors:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f8929b58", + "metadata": {}, + "outputs": [], + "source": [ + "if hasattr(subnet, \"res_cluster_est\"):\n", + " display(\"Estimated load allocation factor = \" + \"{:.2f}\".format(subnet.res_cluster_est[0]))\n", + " display(\"Estimated generation allocation factor = \" + \"{:.2f}\".format(subnet.res_cluster_est[1]))" + ] + }, + { + "cell_type": "markdown", + "id": "f540c85f", + "metadata": {}, + "source": [ + "### State Estimation, Example 2 - LPN grid with realistic measurement uncertainty\n", + "\n", + "In this example, state estimation will be shown on the same portion of the LPN grid used in example 1, but measurement uncertainties will be modified to reflect more realistic conditions usually present in real grids." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "da2533ac", + "metadata": {}, + "outputs": [], + "source": [ + "# Clean measurements already present in the grid\n", + "subnet.measurement.drop(subnet.measurement.index, inplace=True)\n", + "\n", + "# Create the voltage measurements with 0.2% uncertainty\n", + "subnet = create_ext_grid_voltage_measurements(subnet, unc=0.2)\n", + "\n", + "# Create the active and reactive power measurements with 1% uncertainty\n", + "subnet = create_ext_grid_inj_measurements(subnet, unc=1)\n", + "\n", + "# Remove injection values from shunt or ward elements (which would not be seen as power injection measurements)\n", + "remove_shunt_injection_from_meas(subnet,\"shunt\")\n", + "remove_shunt_injection_from_meas(subnet,\"ward\")" + ] + }, + { + "cell_type": "markdown", + "id": "a56f60cd", + "metadata": {}, + "source": [ + "Visualize the measurements used for state estimation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3ecefe18", + "metadata": {}, + "outputs": [], + "source": [ + "display(subnet.measurement)" + ] + }, + { + "cell_type": "markdown", + "id": "098f3566", + "metadata": {}, + "source": [ + "Run state estimation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0b3e799d", + "metadata": {}, + "outputs": [], + "source": [ + "# Run the state estimation algorithm\n", + "try: \n", + " success = se.estimate(subnet, algorithm=\"af-wls\", tolerance=1e-4, maximum_iterations=500)\n", + " display(\"State estimation successfully converged in \" + str(success[\"num_iterations\"]) + \" iterations.\")\n", + "except: \n", + " display(\"State estimation did not converge\")" + ] + }, + { + "cell_type": "markdown", + "id": "909e3380", + "metadata": {}, + "source": [ + "Compare the voltage magnitude state estimation results to the reference results given by the initial power flow; in this case, higher errors should be expected with respect to example 1, due to the larger uncertainty of the measurements. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "05335535", + "metadata": {}, + "outputs": [], + "source": [ + "if subnet.res_bus_est.index.size > 1: \n", + " volt_diff = subnet.res_bus[\"vm_pu\"].values - subnet.res_bus_est[\"vm_pu\"].values\n", + " max_volt_diff = np.max(abs(volt_diff))\n", + " display(\"Maximum voltage difference = \" + \"{:.6f}\".format(max_volt_diff) + \" p.u.\")" + ] + }, + { + "cell_type": "markdown", + "id": "20ecc507", + "metadata": {}, + "source": [ + "Compare line current magnitude state estimation results to reference results given by the initial power flow; in this case, higher errors should be expected with respect to example 1, due to the larger uncertainty of the measurements. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3af6cf54", + "metadata": {}, + "outputs": [], + "source": [ + "if subnet.res_line_est.index.size > 1: \n", + " curr_diff = subnet.res_line[\"i_from_ka\"].values - subnet.res_line_est[\"i_from_ka\"].values\n", + " max_curr_diff = np.max(abs(curr_diff))\n", + " display(\"Maximum current difference = \" + \"{:.6f}\".format(max_curr_diff) + \" kA\")" + ] + }, + { + "cell_type": "markdown", + "id": "52983710", + "metadata": {}, + "source": [ + "Visualize the estimated allocation factors:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8c65f72e", + "metadata": {}, + "outputs": [], + "source": [ + "if hasattr(subnet, \"res_cluster_est\"):\n", + " display(\"Estimated load allocation factor = \" + \"{:.2f}\".format(subnet.res_cluster_est[0]))\n", + " display(\"Estimated generation allocation factor = \" + \"{:.2f}\".format(subnet.res_cluster_est[1]))" + ] + }, + { + "cell_type": "markdown", + "id": "a55ff0c9", + "metadata": {}, + "source": [ + "### State Estimation, Example 3 - SPN grid (realistic measurement uncertainty)\n", + "\n", + "In this example, state estimation will be shown on a portion of the SPN grid. Realistic measurement uncertainties will be considered. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9bef30de", + "metadata": {}, + "outputs": [], + "source": [ + "# Import the grid for the analysis\n", + "filename = \"SPN EQ SSH_0401_eq.json\" # Give here the name of the json file with the UKPN grid you want to use\n", + "if os.path.isfile(filename):\n", + " net = pp.from_json(filename)\n", + "else:\n", + " print(\"file does not exist, creating a dummy net\")\n", + " net = pp.create_empty_network()\n", + " bus = pp.create_bus(net, vn_kv=132)\n", + " pp.create_ext_grid(net, bus=bus)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e46eaa48", + "metadata": {}, + "outputs": [], + "source": [ + "# Apply the workarounds on the selected grid\n", + "license_area = \"SPN\" # Provide here the name of the considered license area. It should be \"LPN\", \"SPN\", or \"EPN\".\n", + "if net.bus.index.size > 1:\n", + " remove_impedance = True # Decide if removing fictious impedances from the grid or not\n", + " net = apply_workarounds(net, license_area, remove_impedance)" + ] + }, + { + "cell_type": "markdown", + "id": "68d1b14e", + "metadata": {}, + "source": [ + "Apply the grid reduction to focus the analysis on a limited portion of the overall network" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ef68d3cc", + "metadata": {}, + "outputs": [], + "source": [ + "start_bus = 4101 # Select the bus of interest around which you want to reduce the grid\n", + "\n", + "# Call the main function for grid reduction\n", + "subnet, trafo_sens_df, trafo3w_sens_df, impedance_sens_df, boundaries, created_ext_grids, created_pq_injections = build_reduced_network(\n", + " net, \n", + " start_bus=start_bus, # start bus considered for the reduction\n", + " method=\"power_flow\", # used method (only power flow available in this tutorial)\n", + " sensitivity_threshold=0.05, # threshold to decide if cutting or not the subnet\n", + " min_working_current_ka=0.001, # minimum current limit considered for the cutting\n", + " vn_max_kv=50.0, # maximum voltage limit considered for the cutting\n", + " deltaP_MW=1.0, # delta of active power toapplied for the sensitivity calculation\n", + " deltaQ_Mvar=0.0, # delta of reactive power toapplied for the sensitivity calculation\n", + " cut_downward_elements=True) # decide if apply cuts also in downstream direction (lower voltage levels) or not\n", + "\n", + "if ~np.any(net.bus.index == start_bus):\n", + " print(\"The selected bus was not found in the considered grid\")\n", + "else:\n", + " print(\"kept buses:\", len(subnet.bus))\n", + " print(\"boundaries:\", len(boundaries[\"hv\"])+len(boundaries[\"lv\"]))\n", + " print(\"ext_grids created:\", len(created_ext_grids))\n", + " print(\"pq_injections_created:\", len(created_pq_injections))" + ] + }, + { + "cell_type": "markdown", + "id": "1d49f1b8", + "metadata": {}, + "source": [ + "Run the power flow to create the conditions of the grid assumed as reference values and extract the measurements from such conditions." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e5adf5f8", + "metadata": {}, + "outputs": [], + "source": [ + "# Run the power flow\n", + "pp.runpp(subnet, run_control=False, lightsim2grid=False, max_iteration=100)\n", + "\n", + "# Clean measurements already present in the grid\n", + "subnet.measurement.drop(subnet.measurement.index, inplace=True)\n", + "\n", + "# Create the voltage measurements with 0.2% uncertainty\n", + "subnet = create_ext_grid_voltage_measurements(subnet, unc=0.2)\n", + "\n", + "# Create the active and reactive power measurements 1% uncertainty\n", + "subnet = create_ext_grid_inj_measurements(subnet, unc=1)\n", + "\n", + "# Remove injection values from shunt or ward elements (which would not be seen as power injection measurements)\n", + "remove_shunt_injection_from_meas(subnet,\"shunt\")\n", + "remove_shunt_injection_from_meas(subnet,\"ward\")" + ] + }, + { + "cell_type": "markdown", + "id": "a3b3fc92", + "metadata": {}, + "source": [ + "Visualize the measurements that will be used for state estimation purposes" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d8f11c77", + "metadata": {}, + "outputs": [], + "source": [ + "display(subnet.measurement)\n", + "display(\"Total nodes in the grid: \" + str(len(subnet.bus)))\n", + "display(\"Total number of measurements: \" + str(len(subnet.measurement)))\n", + "display(\"Measurement redundancy: \" + \"{:.2f}\".format(100*len(subnet.measurement)/(2*len(subnet.bus))) + \" %\")" + ] + }, + { + "cell_type": "markdown", + "id": "974187e0", + "metadata": {}, + "source": [ + "Create the load and generation clusters needed for the AF-WLS algorithm" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0cd94854", + "metadata": {}, + "outputs": [], + "source": [ + "# Create the load and generation clusters under the load and sgen \"type\" attribute\n", + "subnet.load[\"type\"] = \"generic_load\"\n", + "subnet.sgen[\"type\"] = \"generic_sgen\"" + ] + }, + { + "cell_type": "markdown", + "id": "0ccf63df", + "metadata": {}, + "source": [ + "Modify the nominal values of P and Q for loads and sgens" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cd5720f5", + "metadata": {}, + "outputs": [], + "source": [ + "# Modify the nominal values of power\n", + "subnet.load[\"p_mw\"] *= 5\n", + "subnet.load[\"q_mvar\"] *= 5\n", + "subnet.sgen[\"p_mw\"] *= 3\n", + "subnet.sgen[\"q_mvar\"] *= 3" + ] + }, + { + "cell_type": "markdown", + "id": "bc2c2331", + "metadata": {}, + "source": [ + "Run state estimation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6302573d", + "metadata": {}, + "outputs": [], + "source": [ + "# Run the state estimation algorithm\n", + "try: \n", + " success = se.estimate(subnet, algorithm=\"af-wls\", tolerance=1e-4, maximum_iterations=500)\n", + " display(\"State estimation successfully converged in \" + str(success[\"num_iterations\"]) + \" iterations.\")\n", + "except: \n", + " display(\"State estimation did not converge\")" + ] + }, + { + "cell_type": "markdown", + "id": "b9c8f9dc", + "metadata": {}, + "source": [ + "Compare voltage magnitude state estimation results to reference results given by the initial power flow:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ee8d2173", + "metadata": {}, + "outputs": [], + "source": [ + "if subnet.res_bus_est.index.size > 1: \n", + " volt_diff = subnet.res_bus[\"vm_pu\"].values - subnet.res_bus_est[\"vm_pu\"].values\n", + " max_volt_diff = np.max(abs(volt_diff))\n", + " display(\"Maximum voltage difference = \" + \"{:.6f}\".format(max_volt_diff) + \" p.u.\")" + ] + }, + { + "cell_type": "markdown", + "id": "6fdbea62", + "metadata": {}, + "source": [ + "Compare line current magnitude state estimation results to reference results given by the initial power flow:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "733b87b8", + "metadata": {}, + "outputs": [], + "source": [ + "if subnet.res_line_est.index.size > 1: \n", + " curr_diff = subnet.res_line[\"i_from_ka\"].values - subnet.res_line_est[\"i_from_ka\"].values\n", + " max_curr_diff = np.max(abs(curr_diff))\n", + " display(\"Maximum current difference = \" + \"{:.6f}\".format(max_curr_diff) + \" kA\")" + ] + }, + { + "cell_type": "markdown", + "id": "375c4546", + "metadata": {}, + "source": [ + "Visualize the estimated allocation factors:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1ab912d2", + "metadata": {}, + "outputs": [], + "source": [ + "if hasattr(subnet, \"res_cluster_est\"):\n", + " display(\"Estimated load allocation factor = \" + \"{:.2f}\".format(subnet.res_cluster_est[0]))\n", + " display(\"Estimated generation allocation factor = \" + \"{:.2f}\".format(subnet.res_cluster_est[1]))" + ] + }, + { + "cell_type": "markdown", + "id": "5d78f02a", + "metadata": {}, + "source": [ + "### State Estimation, Example 4 - EPN grid (realistic measurement uncertainty)\n", + "\n", + "In this example, state estimation will be shown on a portion of the EPN grid. Realistic measurement uncertainties will be considered. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4c693b53", + "metadata": {}, + "outputs": [], + "source": [ + "# Import the grid for the analysis\n", + "filename = \"EPN EQ SSH_0401_eq.json\" # Give here the name of the json file with the UKPN grid you want to use\n", + "if os.path.isfile(filename):\n", + " net = pp.from_json(filename)\n", + "else:\n", + " print(\"file does not exist, creating a dummy net\")\n", + " net = pp.create_empty_network()\n", + " bus = pp.create_bus(net, vn_kv=132)\n", + " pp.create_ext_grid(net, bus=bus)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ee941f07", + "metadata": {}, + "outputs": [], + "source": [ + "# Apply the workarounds on the selected grid\n", + "license_area = \"EPN\" # Provide here the name of the considered license area. It should be \"LPN\", \"SPN\", or \"EPN\".\n", + "if net.bus.index.size > 1:\n", + " remove_impedance = True # Decide if removing fictious impedances from the grid or not\n", + " net = apply_workarounds(net, license_area, remove_impedance)" + ] + }, + { + "cell_type": "markdown", + "id": "2d609375", + "metadata": {}, + "source": [ + "Apply the grid reduction to focus the analysis on a limited portion of the overall network" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d91d46b9", + "metadata": {}, + "outputs": [], + "source": [ + "start_bus = 150 # Select the bus of interest around which you want to reduce the grid\n", + "\n", + "# Call the main function for grid reduction\n", + "subnet, trafo_sens_df, trafo3w_sens_df, impedance_sens_df, boundaries, created_ext_grids, created_pq_injections = build_reduced_network(\n", + " net, \n", + " start_bus=start_bus, # start bus considered for the reduction\n", + " method=\"power_flow\", # used method (only power flow available in this tutorial)\n", + " sensitivity_threshold=0.05, # threshold to decide if cutting or not the subnet\n", + " min_working_current_ka=0.001, # minimum current limit considered for the cutting\n", + " vn_max_kv=50.0, # maximum voltage limit considered for the cutting\n", + " deltaP_MW=1.0, # delta of active power toapplied for the sensitivity calculation\n", + " deltaQ_Mvar=0.0, # delta of reactive power toapplied for the sensitivity calculation\n", + " cut_downward_elements=True) # decide if apply cuts also in downstream direction (lower voltage levels) or not\n", + "\n", + "if ~np.any(net.bus.index == start_bus):\n", + " print(\"The selected bus was not found in the considered grid\")\n", + "else:\n", + " print(\"kept buses:\", len(subnet.bus))\n", + " print(\"boundaries:\", len(boundaries[\"hv\"])+len(boundaries[\"lv\"]))\n", + " print(\"ext_grids created:\", len(created_ext_grids))\n", + " print(\"pq_injections_created:\", len(created_pq_injections))" + ] + }, + { + "cell_type": "markdown", + "id": "57fb7d37", + "metadata": {}, + "source": [ + "Run the power flow to create the conditions of the grid assumed as reference values and extract the measurements from such conditions." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e1cb085b", + "metadata": {}, + "outputs": [], + "source": [ + "# Run the power flow\n", + "pp.runpp(subnet, run_control=False, lightsim2grid=False, max_iteration=100)\n", + "\n", + "# Clean measurements already present in the grid\n", + "subnet.measurement.drop(subnet.measurement.index, inplace=True)\n", + "\n", + "# Create the voltage measurements with 0.2% uncertainty\n", + "subnet = create_ext_grid_voltage_measurements(subnet, unc=0.2)\n", + "\n", + "# Create the active and reactive power measurements with 1% uncertainty\n", + "subnet = create_ext_grid_inj_measurements(subnet, unc=1)\n", + "\n", + "# Remove injection values from shunt or ward elements (which would not be seen as power injection measurements)\n", + "remove_shunt_injection_from_meas(subnet,\"shunt\")\n", + "remove_shunt_injection_from_meas(subnet,\"ward\")" + ] + }, + { + "cell_type": "markdown", + "id": "9df1104b", + "metadata": {}, + "source": [ + "Visualize the measurements that will be used for state estimation purposes" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "031c7c4d", + "metadata": {}, + "outputs": [], + "source": [ + "display(subnet.measurement)\n", + "display(\"Total nodes in the grid: \" + str(len(subnet.bus)))\n", + "display(\"Total number of measurements: \" + str(len(subnet.measurement)))\n", + "display(\"Measurement redundancy: \" + \"{:.2f}\".format(100*len(subnet.measurement)/(2*len(subnet.bus))) + \" %\")" + ] + }, + { + "cell_type": "markdown", + "id": "0e5f1f58", + "metadata": {}, + "source": [ + "Create the load and generation clusters needed for the AF-WLS algorithm" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "137ed26f", + "metadata": {}, + "outputs": [], + "source": [ + "# Create the load and generation clusters under the load and sgen \"type\" attribute\n", + "subnet.load[\"type\"] = \"generic_load\"\n", + "subnet.sgen[\"type\"] = \"generic_sgen\"" + ] + }, + { + "cell_type": "markdown", + "id": "d02f7526", + "metadata": {}, + "source": [ + "Modify the nominal values of P and Q for loads and sgens" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e221a8e3", + "metadata": {}, + "outputs": [], + "source": [ + "# Modify the nominal values of power\n", + "subnet.load[\"p_mw\"] *= 1.5\n", + "subnet.load[\"q_mvar\"] *= 1.5\n", + "subnet.sgen[\"p_mw\"] *= 1.2\n", + "subnet.sgen[\"q_mvar\"] *= 1.2" + ] + }, + { + "cell_type": "markdown", + "id": "00ec41b3", + "metadata": {}, + "source": [ + "Run state estimation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e4ef2fd2", + "metadata": {}, + "outputs": [], + "source": [ + "# Run the state estimation algorithm\n", + "try: \n", + " success = se.estimate(subnet, algorithm=\"af-wls\", tolerance=1e-4, maximum_iterations=500)\n", + " display(\"State estimation successfully converged in \" + str(success[\"num_iterations\"]) + \" iterations.\")\n", + "except: \n", + " display(\"State estimation did not converge\")" + ] + }, + { + "cell_type": "markdown", + "id": "30e29e4c", + "metadata": {}, + "source": [ + "Compare voltage magnitude state estimation results to reference results given by the initial power flow:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cd57ca9a", + "metadata": {}, + "outputs": [], + "source": [ + "if subnet.res_bus_est.index.size > 1: \n", + " volt_diff = subnet.res_bus[\"vm_pu\"].values - subnet.res_bus_est[\"vm_pu\"].values\n", + " max_volt_diff = np.max(abs(volt_diff))\n", + " display(\"Maximum voltage difference = \" + \"{:.6f}\".format(max_volt_diff) + \" p.u.\")" + ] + }, + { + "cell_type": "markdown", + "id": "790895af", + "metadata": {}, + "source": [ + "Compare line current magnitude state estimation results to reference results given by the initial power flow:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "512898fa", + "metadata": {}, + "outputs": [], + "source": [ + "if subnet.res_line_est.index.size > 1: \n", + " curr_diff = subnet.res_line[\"i_from_ka\"].values - subnet.res_line_est[\"i_from_ka\"].values\n", + " max_curr_diff = np.max(abs(curr_diff))\n", + " display(\"Maximum current difference = \" + \"{:.6f}\".format(max_curr_diff) + \" kA\")" + ] + }, + { + "cell_type": "markdown", + "id": "e9bf7223", + "metadata": {}, + "source": [ + "Visualize the estimated allocation factors:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "70600317", + "metadata": {}, + "outputs": [], + "source": [ + "if hasattr(subnet, \"res_cluster_est\"):\n", + " display(\"Estimated load allocation factor = \" + \"{:.2f}\".format(subnet.res_cluster_est[0]))\n", + " display(\"Estimated generation allocation factor = \" + \"{:.2f}\".format(subnet.res_cluster_est[1]))" + ] + }, + { + "cell_type": "markdown", + "id": "461d3dd0", + "metadata": {}, + "source": [ + "### State Forecasting, Example 1 - LPN Grid\n", + "\n", + "State forecasting differs from state estimation because forecast data, instead of real-time measurements, are used to estimate the future operating conditions of the grid. \n", + "The way state forecasting is performed strictly depends on the available forecast data. \n", + "\n", + "A first option, it to use the time series of real-time measurements for predicting future measurements (same values that are measured - and at the same location - but in the future).\n", + "In this case, the same AF-WLS approach seen in the previous examples for state estimation can be used (the only difference is that the values present in the net.measurement table will be forecast values instead of real measurements).\n", + "\n", + "Since forecast data can be generated in multiple ways, also from other data rather than measurement time series, grid operators often have the availability of predicted consumption or generation at the different buses of the grid. \n", + "In this case, a classical Weighted Least Squares algorithm can be used if the forecast data are enough to reach the full observability of the grid. \n", + "The following examples will show the use of the WLS formulation for state forecasting in the case of fully observable grids. \n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2ef93a9c", + "metadata": {}, + "outputs": [], + "source": [ + "# Import the grid for the analysis\n", + "filename = \"LPN EQ SSH_0401_eq.json\" # Give here the name of the json file with the UKPN grid you want to use\n", + "if os.path.isfile(filename):\n", + " net = pp.from_json(filename)\n", + "else:\n", + " print(\"file does not exist, creating a dummy net\")\n", + " net = pp.create_empty_network()\n", + " bus = pp.create_bus(net, vn_kv=132)\n", + " pp.create_ext_grid(net, bus=bus)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4433ccdf", + "metadata": {}, + "outputs": [], + "source": [ + "# Apply the workarounds on the selected grid\n", + "license_area = \"LPN\" # Provide here the name of the considered license area. It should be \"LPN\", \"SPN\", or \"EPN\".\n", + "if net.bus.index.size > 1:\n", + " remove_impedance = True # Decide if removing fictious impedances from the grid or not\n", + " net = apply_workarounds(net, license_area, remove_impedance)" + ] + }, + { + "cell_type": "markdown", + "id": "b9cdfb4b", + "metadata": {}, + "source": [ + "Apply the grid reduction to focus the analysis on a limited portion of the overall grid" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "14d78cef", + "metadata": {}, + "outputs": [], + "source": [ + "start_bus = 1712 # Select the bus of interest around which you want to reduce the grid\n", + "\n", + "# Call the main function for grid reduction\n", + "subnet, trafo_sens_df, trafo3w_sens_df, impedance_sens_df, boundaries, created_ext_grids, created_pq_injections = build_reduced_network(\n", + " net, \n", + " start_bus=start_bus, # start bus considered for the reduction\n", + " method=\"power_flow\", # used method (only power flow available in this tutorial)\n", + " sensitivity_threshold=0.05, # threshold to decide if cutting or not the subnet\n", + " min_working_current_ka=0.001, # minimum current limit considered for the cutting\n", + " vn_max_kv=50.0, # maximum voltage limit considered for the cutting\n", + " deltaP_MW=1.0, # delta of active power toapplied for the sensitivity calculation\n", + " deltaQ_Mvar=0.0, # delta of reactive power toapplied for the sensitivity calculation\n", + " cut_downward_elements=True) # decide if apply cuts also in downstream direction (lower voltage levels) or not\n", + "\n", + "if ~np.any(net.bus.index == start_bus):\n", + " print(\"The selected bus was not found in the considered grid\")\n", + "else:\n", + " print(\"kept buses:\", len(subnet.bus))\n", + " print(\"boundaries:\", len(boundaries[\"hv\"])+len(boundaries[\"lv\"]))\n", + " print(\"ext_grids created:\", len(created_ext_grids))\n", + " print(\"pq_injections_created:\", len(created_pq_injections))" + ] + }, + { + "cell_type": "markdown", + "id": "059e6a70", + "metadata": {}, + "source": [ + "Run a power flow to create the conditions of the grid assumed as reference values and extract forecasted values (with uncertainty) from such conditions." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c48d2a30", + "metadata": {}, + "outputs": [], + "source": [ + "# Run the power flow\n", + "pp.runpp(subnet, run_control=False, lightsim2grid=False, max_iteration=100)\n", + "\n", + "# Clean measurements already present in the grid\n", + "subnet.measurement.drop(subnet.measurement.index, inplace=True)\n", + "\n", + "# Create the voltage measurements with 0.5% uncertainty\n", + "subnet = create_ext_grid_voltage_measurements(subnet, unc=0.5)\n", + "\n", + "# Create the active and reactive power measurements with 10% uncertainty\n", + "subnet = create_bus_inj_measurements(subnet, unc=10)\n", + "\n", + "# Remove injection values from shunt or ward elements (which would not be seen as power injection measurements)\n", + "remove_shunt_injection_from_meas(subnet,\"shunt\")\n", + "remove_shunt_injection_from_meas(subnet,\"ward\")" + ] + }, + { + "cell_type": "markdown", + "id": "73510ee5", + "metadata": {}, + "source": [ + "Visualize the data that will be used for state forecasting" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4a3c5beb", + "metadata": {}, + "outputs": [], + "source": [ + "display(\"Total nodes in the grid: \" + str(len(subnet.bus)))\n", + "display(\"Total number of forecasts: \" + str(len(subnet.measurement)))\n", + "display(\"Data redundancy: \" + \"{:.2f}\".format(100*len(subnet.measurement)/(2*len(subnet.bus))) + \" %\")" + ] + }, + { + "cell_type": "markdown", + "id": "e11f43da", + "metadata": {}, + "source": [ + "Run the WLS algorithm" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5f947c3e", + "metadata": {}, + "outputs": [], + "source": [ + "# Run the state forecasting algorithm\n", + "try: \n", + " success = se.estimate(subnet, algorithm=\"wls\", tolerance=1e-6, maximum_iterations=100)\n", + " display(\"State estimation successfully converged in \" + str(success[\"num_iterations\"]) + \" iterations.\")\n", + "except: \n", + " display(\"State estimation did not converge\")" + ] + }, + { + "cell_type": "markdown", + "id": "47cba180", + "metadata": {}, + "source": [ + "Compare voltage magnitude state forecasting results to the reference values given by the initial power flow:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9359f0d4", + "metadata": {}, + "outputs": [], + "source": [ + "if net.res_bus_est.index.size > 1: \n", + " volt_diff = subnet.res_bus[\"vm_pu\"].values - subnet.res_bus_est[\"vm_pu\"].values\n", + " max_volt_diff = np.max(abs(volt_diff))\n", + " display(\"Maximum voltage difference = \" + \"{:.6f}\".format(max_volt_diff) + \" p.u.\")" + ] + }, + { + "cell_type": "markdown", + "id": "a998a20b", + "metadata": {}, + "source": [ + "Compare line current magnitude state forecasting results to the reference values given by the power flow" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fb51c78b", + "metadata": {}, + "outputs": [], + "source": [ + "if subnet.res_line_est.index.size > 1: \n", + " curr_diff = subnet.res_line[\"i_from_ka\"].values - subnet.res_line_est[\"i_from_ka\"].values\n", + " max_curr_diff = np.max(abs(curr_diff))\n", + " display(\"Maximum current difference = \" + \"{:.6f}\".format(max_curr_diff) + \" kA\")" + ] + }, + { + "cell_type": "markdown", + "id": "6484f361", + "metadata": {}, + "source": [ + "### State Forecasting, Example 2 - SPN Grid\n", + "\n", + "In this example, state forecasting will be shown on a portion of the SPN grid." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a95824b1", + "metadata": {}, + "outputs": [], + "source": [ + "# Import the grid for the analysis\n", + "filename = \"SPN EQ SSH_0401_eq.json\" # Give here the name of the json file with the UKPN grid you want to use\n", + "if os.path.isfile(filename):\n", + " net = pp.from_json(filename)\n", + "else:\n", + " print(\"file does not exist, creating a dummy net\")\n", + " net = pp.create_empty_network()\n", + " bus = pp.create_bus(net, vn_kv=132)\n", + " pp.create_ext_grid(net, bus=bus)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "445c0b7e", + "metadata": {}, + "outputs": [], + "source": [ + "# Apply the workarounds on the selected grid\n", + "license_area = \"SPN\" # Provide here the name of the considered license area. It should be \"LPN\", \"SPN\", or \"EPN\".\n", + "if net.bus.index.size > 1:\n", + " remove_impedance = True # Decide if removing fictious impedances from the grid or not\n", + " net = apply_workarounds(net, license_area, remove_impedance)" + ] + }, + { + "cell_type": "markdown", + "id": "c42653cd", + "metadata": {}, + "source": [ + "Apply the grid reduction to focus the analysis on a limited portion of the overall grid" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fd25e3ff", + "metadata": {}, + "outputs": [], + "source": [ + "start_bus = 4101 # Select the bus of interest around which you want to reduce the grid\n", + "\n", + "# Call the main function for grid reduction\n", + "subnet, trafo_sens_df, trafo3w_sens_df, impedance_sens_df, boundaries, created_ext_grids, created_pq_injections = build_reduced_network(\n", + " net, \n", + " start_bus=start_bus, # start bus considered for the reduction\n", + " method=\"power_flow\", # used method (only power flow available in this tutorial)\n", + " sensitivity_threshold=0.05, # threshold to decide if cutting or not the subnet\n", + " min_working_current_ka=0.001, # minimum current limit considered for the cutting\n", + " vn_max_kv=50.0, # maximum voltage limit considered for the cutting\n", + " deltaP_MW=1.0, # delta of active power toapplied for the sensitivity calculation\n", + " deltaQ_Mvar=0.0, # delta of reactive power toapplied for the sensitivity calculation\n", + " cut_downward_elements=True) # decide if apply cuts also in downstream direction (lower voltage levels) or not\n", + "\n", + "if ~np.any(net.bus.index == start_bus):\n", + " print(\"The selected bus was not found in the considered grid\")\n", + "else:\n", + " print(\"kept buses:\", len(subnet.bus))\n", + " print(\"boundaries:\", len(boundaries[\"hv\"])+len(boundaries[\"lv\"]))\n", + " print(\"ext_grids created:\", len(created_ext_grids))\n", + " print(\"pq_injections_created:\", len(created_pq_injections))" + ] + }, + { + "cell_type": "markdown", + "id": "37abc657", + "metadata": {}, + "source": [ + "Run a power flow to create the conditions of the grid assumed as reference values and extract forecasted values (with uncertainty) from such conditions." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ffde83cd", + "metadata": {}, + "outputs": [], + "source": [ + "# Run the power flow\n", + "pp.runpp(subnet, run_control=False, lightsim2grid=False, max_iteration=100)\n", + "\n", + "# Clean measurements already present in the grid\n", + "subnet.measurement.drop(subnet.measurement.index, inplace=True)\n", + "\n", + "# Create the voltage measurements with 0.5% uncertainty\n", + "subnet = create_ext_grid_voltage_measurements(subnet, unc=0.5)\n", + "\n", + "# Create the active and reactive power measurements with 10% uncertainty\n", + "subnet = create_bus_inj_measurements(subnet, unc=10)\n", + "\n", + "# Remove injection values from shunt or ward elements (which would not be seen as power injection measurements)\n", + "remove_shunt_injection_from_meas(subnet,\"shunt\")\n", + "remove_shunt_injection_from_meas(subnet,\"ward\")" + ] + }, + { + "cell_type": "markdown", + "id": "8fdfa103", + "metadata": {}, + "source": [ + "Visualize the data that will be used for state forecasting" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "099766b5", + "metadata": {}, + "outputs": [], + "source": [ + "display(\"Total nodes in the grid: \" + str(len(subnet.bus)))\n", + "display(\"Total number of forecasts: \" + str(len(subnet.measurement)))\n", + "display(\"Data redundancy: \" + \"{:.2f}\".format(100*len(subnet.measurement)/(2*len(subnet.bus))) + \" %\")" + ] + }, + { + "cell_type": "markdown", + "id": "4a638a22", + "metadata": {}, + "source": [ + "Run the WLS algorithm" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c270643d", + "metadata": {}, + "outputs": [], + "source": [ + "# Run the state forecasting algorithm\n", + "try: \n", + " success = se.estimate(subnet, algorithm=\"wls\", tolerance=1e-6, maximum_iterations=100)\n", + " display(\"State estimation successfully converged in \" + str(success[\"num_iterations\"]) + \" iterations.\")\n", + "except: \n", + " display(\"State estimation did not converge\")" + ] + }, + { + "cell_type": "markdown", + "id": "72d84c1f", + "metadata": {}, + "source": [ + "Compare voltage magnitude state forecasting results to the reference values given by the initial power flow:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0a03fc28", + "metadata": {}, + "outputs": [], + "source": [ + "if subnet.res_bus_est.index.size > 1: \n", + " volt_diff = subnet.res_bus[\"vm_pu\"].values - subnet.res_bus_est[\"vm_pu\"].values\n", + " max_volt_diff = np.max(abs(volt_diff))\n", + " display(\"Maximum voltage difference = \" + \"{:.6f}\".format(max_volt_diff) + \" p.u.\")" + ] + }, + { + "cell_type": "markdown", + "id": "f7802877", + "metadata": {}, + "source": [ + "Compare line current magnitude state forecasting results to the reference values given by the power flow" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8078fdf6", + "metadata": {}, + "outputs": [], + "source": [ + "if subnet.res_line_est.index.size > 1: \n", + " curr_diff = subnet.res_line[\"i_from_ka\"].values - subnet.res_line_est[\"i_from_ka\"].values\n", + " max_curr_diff = np.max(abs(curr_diff))\n", + " display(\"Maximum current difference = \" + \"{:.6f}\".format(max_curr_diff) + \" kA\")" + ] + }, + { + "cell_type": "markdown", + "id": "1eaa1e75", + "metadata": {}, + "source": [ + "### State Forecasting, Example 3 - EPN Grid\n", + "\n", + "In this example, state forecasting will be shown on a portion of the EPN grid." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f4ef3cc2", + "metadata": {}, + "outputs": [], + "source": [ + "# Import the grid for the analysis\n", + "filename = \"EPN EQ SSH_0401_eq.json\" # Give here the name of the json file with the UKPN grid you want to use\n", + "if os.path.isfile(filename):\n", + " net = pp.from_json(filename)\n", + "else:\n", + " print(\"file does not exist, creating a dummy net\")\n", + " net = pp.create_empty_network()\n", + " bus = pp.create_bus(net, vn_kv=132)\n", + " pp.create_ext_grid(net, bus=bus)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1956f03a", + "metadata": {}, + "outputs": [], + "source": [ + "# Apply the workarounds on the selected grid\n", + "license_area = \"EPN\" # Provide here the name of the considered license area. It should be \"LPN\", \"SPN\", or \"EPN\".\n", + "if net.bus.index.size > 1:\n", + " remove_impedance = True # Decide if removing fictious impedances from the grid or not\n", + " net = apply_workarounds(net, license_area, remove_impedance)" + ] + }, + { + "cell_type": "markdown", + "id": "d56c9cad", + "metadata": {}, + "source": [ + "Apply the grid reduction to focus the analysis on a limited portion of the overall grid" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8b030220", + "metadata": {}, + "outputs": [], + "source": [ + "start_bus = 150 # Select the bus of interest around which you want to reduce the grid\n", + "\n", + "# Call the main function for grid reduction\n", + "subnet, trafo_sens_df, trafo3w_sens_df, impedance_sens_df, boundaries, created_ext_grids, created_pq_injections = build_reduced_network(\n", + " net, \n", + " start_bus=start_bus, # start bus considered for the reduction\n", + " method=\"power_flow\", # used method (only power flow available in this tutorial)\n", + " sensitivity_threshold=0.05, # threshold to decide if cutting or not the subnet\n", + " min_working_current_ka=0.001, # minimum current limit considered for the cutting\n", + " vn_max_kv=50.0, # maximum voltage limit considered for the cutting\n", + " deltaP_MW=1.0, # delta of active power toapplied for the sensitivity calculation\n", + " deltaQ_Mvar=0.0, # delta of reactive power toapplied for the sensitivity calculation\n", + " cut_downward_elements=True) # decide if apply cuts also in downstream direction (lower voltage levels) or not\n", + "\n", + "if ~np.any(net.bus.index == start_bus):\n", + " print(\"The selected bus was not found in the considered grid\")\n", + "else:\n", + " print(\"kept buses:\", len(subnet.bus))\n", + " print(\"boundaries:\", len(boundaries[\"hv\"])+len(boundaries[\"lv\"]))\n", + " print(\"ext_grids created:\", len(created_ext_grids))\n", + " print(\"pq_injections_created:\", len(created_pq_injections))" + ] + }, + { + "cell_type": "markdown", + "id": "5a5e0ef7", + "metadata": {}, + "source": [ + "Run a power flow to create the conditions of the grid assumed as reference values and extract forecasted values (with uncertainty) from such conditions." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "37addf9f", + "metadata": {}, + "outputs": [], + "source": [ + "# Run the power flow\n", + "pp.runpp(subnet, run_control=False, lightsim2grid=False, max_iteration=100)\n", + "\n", + "# Clean measurements already present in the grid\n", + "subnet.measurement.drop(subnet.measurement.index, inplace=True)\n", + "\n", + "# Create the voltage measurements with 0.5% uncertainty\n", + "subnet = create_ext_grid_voltage_measurements(subnet, unc=0.5)\n", + "\n", + "# Create the active and reactive power measurements with 10% uncertainty\n", + "subnet = create_bus_inj_measurements(subnet, unc=10)\n", + "\n", + "# Remove injection values from shunt or ward elements (which would not be seen as power injection measurements)\n", + "remove_shunt_injection_from_meas(subnet,\"shunt\")\n", + "remove_shunt_injection_from_meas(subnet,\"ward\")" + ] + }, + { + "cell_type": "markdown", + "id": "a296ceeb", + "metadata": {}, + "source": [ + "Visualize the data that will be used for state forecasting" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f61bf529", + "metadata": {}, + "outputs": [], + "source": [ + "display(\"Total nodes in the grid: \" + str(len(subnet.bus)))\n", + "display(\"Total number of forecasts: \" + str(len(subnet.measurement)))\n", + "display(\"Data redundancy: \" + \"{:.2f}\".format(100*len(subnet.measurement)/(2*len(subnet.bus))) + \" %\")" + ] + }, + { + "cell_type": "markdown", + "id": "6077e8e1", + "metadata": {}, + "source": [ + "Run the WLS algorithm" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d2014110", + "metadata": {}, + "outputs": [], + "source": [ + "# Run the state forecasting algorithm\n", + "try: \n", + " success = se.estimate(subnet, algorithm=\"wls\", tolerance=1e-6, maximum_iterations=100)\n", + " display(\"State estimation successfully converged in \" + str(success[\"num_iterations\"]) + \" iterations.\")\n", + "except: \n", + " display(\"State estimation did not converge\")" + ] + }, + { + "cell_type": "markdown", + "id": "dd4bd65a", + "metadata": {}, + "source": [ + "Compare voltage magnitude state forecasting results to the reference values given by the initial power flow:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "686d0835", + "metadata": {}, + "outputs": [], + "source": [ + "if subnet.res_bus_est.index.size > 1: \n", + " volt_diff = subnet.res_bus[\"vm_pu\"].values - subnet.res_bus_est[\"vm_pu\"].values\n", + " max_volt_diff = np.max(abs(volt_diff))\n", + " display(\"Maximum voltage difference = \" + \"{:.6f}\".format(max_volt_diff) + \" p.u.\")" + ] + }, + { + "cell_type": "markdown", + "id": "f08b6f0f", + "metadata": {}, + "source": [ + "Compare line current magnitude state forecasting results to the reference values given by the power flow" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "eaec30d6", + "metadata": {}, + "outputs": [], + "source": [ + "if subnet.res_line_est.index.size > 1: \n", + " curr_diff = subnet.res_line[\"i_from_ka\"].values - subnet.res_line_est[\"i_from_ka\"].values\n", + " max_curr_diff = np.max(abs(curr_diff))\n", + " display(\"Maximum current difference = \" + \"{:.6f}\".format(max_curr_diff) + \" kA\")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "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.10.3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}