diff --git a/notebooks/network_analysis_myopic.ipynb b/notebooks/network_analysis_myopic.ipynb new file mode 100644 index 0000000..df61d39 --- /dev/null +++ b/notebooks/network_analysis_myopic.ipynb @@ -0,0 +1,516 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Compare the use of Myopic optimization from Pypsa-Earth-Sec\n", + "\n", + "This notebook aims at exploring the use of Myopic optimization from Pypsa-Earth-Sec.\n", + "The notebook reads:\n", + "- results of myopic optimization for a given list of planning horizons\n", + "- the corresponding overnight optimization results for the same scenarios\n", + "\n", + "Appropriate plots are provided to show at each planning horizon:\n", + "- installed capacity, thus the capacity of power plants already built\n", + "- optimal capacity, thus the total optimized capacity of power plants\n", + "- capacity expansion, which is the difference between optimal and installed capacity\n", + "- energy balance\n", + "- the retired capacity at each time step" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Sources: \n", + "- network_analysis.ipynb: https://github.com/pypsa-meets-earth/documentation/blob/main/notebooks/viz/network_analysis.ipynb\n", + "- Plot capacity - map view: https://github.com/pypsa-meets-earth/documentation/blob/main/notebooks/viz/regional_transm_system_viz.ipynb\n", + "- Analyse energy potential: https://github.com/pypsa-meets-earth/documentation/blob/main/notebooks/build_renewable_profiles.ipynb\n", + "- Analyse energy generation: https://pypsa.readthedocs.io/en/latest/examples/statistics.html" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The execution of the overnight optimizations are assumed to be performed using scenario management with the format, with the format `{country}_{planning_horizon}`.\n", + "The myopic and overnight optimizations are assumed to be performed into two different main repository, which optionally can be the same folder." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "root_myopic = \"/data/davidef/gitsec_my/pypsa-earth-sec/\"\n", + "root_overnight = \"/data/davidef/gitsec_base/pypsa-earth-sec/\"" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Import packages" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import yaml\n", + "import pypsa\n", + "import warnings\n", + "import matplotlib.pyplot as plt\n", + "import geopandas as gpd\n", + "import numpy as np\n", + "import pandas as pd\n", + "from pathlib import Path\n", + "import seaborn as sns\n", + "from datetime import datetime\n", + "from cartopy import crs as ccrs\n", + "from pypsa.plot import add_legend_circles, add_legend_lines, add_legend_patches\n", + "import os\n", + "import xarray as xr\n", + "import cartopy" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# change current directory to parent folder\n", + "import os\n", + "import sys\n", + "\n", + "if not os.path.isdir(\"pypsa-earth\"):\n", + " os.chdir(\"../..\")\n", + "sys.path.append(os.getcwd()+\"/pypsa-earth/scripts\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Path settings\n", + "This section reads the config parameters from your config.yaml file and automatically reads the output of the optimization with those settings" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "config = yaml.safe_load(open(root_myopic + \"config.yaml\"))\n", + "\n", + "# Read config.yaml settings:\n", + "name_myopic = config[\"run\"]\n", + "countries = config[\"countries\"]\n", + "simpl = config[\"scenario\"][\"simpl\"]\n", + "clusters = config[\"scenario\"][\"clusters\"]\n", + "planning_horizons = config[\"scenario\"][\"planning_horizons\"]\n", + "ll = config[\"scenario\"][\"ll\"]\n", + "opts = config[\"scenario\"][\"opts\"]\n", + "sopts = config[\"scenario\"][\"sopts\"]\n", + "demand = config[\"scenario\"][\"demand\"]\n", + "h2export = config[\"export\"][\"h2export\"]\n", + "discountrate = config[\"costs\"][\"discountrate\"]\n", + "\n", + "# Ensure elements are strings and properly joined\n", + "simpl_str = str(simpl[0])\n", + "clusters_str = str(clusters[0])\n", + "ll_str = str(ll[0])\n", + "opts_str = str(opts[0])\n", + "sopts_str = str(sopts[0])\n", + "demand_str = str(demand[0])\n", + "h2export_str = str(h2export[0])\n", + "discountrate_str = str(discountrate[0])\n", + "\n", + "# auxiliary function to get the scenario name using the myopic one as standard.\n", + "# For example: for planning_horizon=2030, NG_MYOPIC becomes NG_2030\n", + "def get_scenario_name(name_myopic, planning_horizon):\n", + " \"\"\"Scenario format {name}_{planning_horizon}\"\"\"\n", + " str_split = name_myopic.split(\"_\")\n", + " return str_split[0] + \"_\" + str(planning_horizon)\n", + "\n", + "# Auxiliary function to get the network path for the myopic configuration\n", + "def get_myopic_network(planning_horizon):\n", + " return os.path.join(\n", + " root_myopic,\n", + " \"results\",\n", + " name_myopic,\n", + " \"postnetworks\",\n", + " f\"elec_s{simpl_str}_{clusters_str}_ec_l{ll_str}_{opts_str}_{sopts_str}_{planning_horizon}_{discountrate_str}_{demand_str}_{h2export_str}export.nc\",\n", + " )\n", + "\n", + "# Auxiliary function to get the network path for the overnight configuration by planning horizon\n", + "def get_overnight_network(planning_horizon):\n", + " name_scenario = get_scenario_name(name_myopic, planning_horizon)\n", + " return os.path.join(\n", + " root_overnight,\n", + " \"results\",\n", + " name_scenario,\n", + " \"postnetworks\", # {clusters_str} , {ll_str} {sopts_str} {h2export_str}\n", + " f\"elec_s{simpl_str}_10_ec_lc1.0_{opts_str}_1H_{planning_horizon}_{discountrate_str}_{demand_str}_10export.nc\",\n", + " )\n", + "\n", + "scenario_name = name_myopic\n", + "scenario_subpath = scenario_name + \"/\" if scenario_name else \"\"\n", + "\n", + "# Country shape file\n", + "regions_onshore_path = root_myopic + f\"pypsa-earth/resources/shapes/country_shapes.geojson\"" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Energy system analysis setup - power and energy generation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "warnings.simplefilter(action='ignore', category=FutureWarning)\n", + "regions_onshore = gpd.read_file(regions_onshore_path)\n", + "country_coordinates = regions_onshore.total_bounds[[0, 2, 1, 3]]\n", + "warnings.simplefilter(action='default', category=FutureWarning)\n", + "\n", + "fp_myopic = get_myopic_network(planning_horizons[0])\n", + "\n", + "n = pypsa.Network(fp_myopic)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Data import check" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Plot of the region of interest" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "x_mean = (regions_onshore.total_bounds[0] + regions_onshore.total_bounds[2])/2\n", + "fig, ax = plt.subplots(figsize=(8, 8), subplot_kw={\"projection\": ccrs.EqualEarth(x_mean)})\n", + "with plt.rc_context({\"patch.linewidth\": 0.}):\n", + " regions_onshore.plot(\n", + " ax=ax,\n", + " facecolor=\"green\",\n", + " edgecolor=\"white\",\n", + " aspect=\"equal\",\n", + " transform=ccrs.PlateCarree(),\n", + " linewidth=0,\n", + " )\n", + "ax.set_title(\", \".join(regions_onshore.name.values))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Analyse energy system" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Analys the future generation capacity expansion of the energy system - bar chart" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Auxiliary function to calculate optimal, installed, and capacity expansion capacities and energy balance\n", + "# for a given network n.\n", + "# It returns a pandas series for each whose name can contain an optional postfix.\n", + "def get_expansion_quantities(n, postfix=\"\"):\n", + " optimal_capacity = n.statistics.optimal_capacity(comps=[\"Generator\"]).droplevel(0).div(1e3).rename(\"optimal\" + postfix).drop(\"load\", errors=\"ignore\")\n", + " installed_capacity = n.statistics.installed_capacity(comps=[\"Generator\"]).droplevel(0).div(1e3).rename(\"installed\" + postfix).drop(\"load\", errors=\"ignore\")\n", + " generation_capacity_expansion = (optimal_capacity - installed_capacity).rename(\"expansion\" + postfix).drop(\"load\", errors=\"ignore\")\n", + " rename_cols = {\n", + " '-': 'Load',\n", + " 'load': 'load shedding',\n", + " }\n", + "\n", + " energy_balance = (\n", + " n.statistics.energy_balance(comps=[\"Generator\", \"Store\"])\n", + " .loc[:, :, \"AC\"]\n", + " .groupby(\"carrier\")\n", + " .sum()\n", + " .div(1e6)\n", + " .rename(\"balance\" + postfix)\n", + " )\n", + " return optimal_capacity, installed_capacity, generation_capacity_expansion, energy_balance" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Loop over the planning horizons and compare the myopic and overnight optimization results" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "scenario_configs = (\n", + " [\"myopic - \" + str(ph) for ph in planning_horizons] +\n", + " [\"overnight - \" + str(ph) for ph in planning_horizons]\n", + ")\n", + "\n", + "optimal_list = []\n", + "installed_list = []\n", + "expansion_list = []\n", + "balance_list = []\n", + "\n", + "for scenario_name in scenario_configs:\n", + "\n", + " config_name, planning_horizon = scenario_name.split(\" - \")\n", + "\n", + " if config_name.lower() == \"myopic\":\n", + " n = pypsa.Network(get_myopic_network(planning_horizon))\n", + " else:\n", + " n = pypsa.Network(get_overnight_network(planning_horizon))\n", + " optimal_capacity, installed_capacity, generation_capacity_expansion, energy_balance = get_expansion_quantities(n, f\" - {scenario_name}\")\n", + " optimal_list.append(optimal_capacity)\n", + " installed_list.append(installed_capacity)\n", + " expansion_list.append(generation_capacity_expansion)\n", + " balance_list.append(energy_balance)\n", + "\n", + "# Merge the pandas series into dataframes\n", + "total_optimal = pd.concat(optimal_list, axis=1)\n", + "total_installed = pd.concat(installed_list, axis=1)\n", + "total_expansion = pd.concat(expansion_list, axis=1)\n", + "total_balance = pd.concat(balance_list, axis=1)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Define an auxiliary function to help the plotting of the dataframe calculated previously about capacity expansion" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def plot_expansion_dataframe(df, title_name):\n", + " plot_df = df.copy().T\n", + "\n", + " colors = {key.lower(): value.lower() for key, value in config[\"plotting\"][\"tech_colors\"].items()}\n", + " colors[\"Combined-Cycle Gas\"] = colors[\"ocgt\"]\n", + " colors[\"load\"] = \"pink\"\n", + "\n", + " # color-matching\n", + " color_list = []\n", + " for col in plot_df.columns:\n", + " original_name = col.lower()\n", + " color = colors.get(original_name.lower(), 'red')\n", + " color_list.append(color)\n", + "\n", + " fig, ax = plt.subplots()\n", + " plot_df.plot.bar(stacked=True, ax=ax, title=title_name, color=color_list)\n", + " handles, labels = ax.get_legend_handles_labels()\n", + " nice_labels = plot_df.columns\n", + " ax.legend(handles, nice_labels, bbox_to_anchor=(1, 0), loc=\"lower left\", title=None, ncol=1)\n", + "\n", + " plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Plot of Expansion of generation capacity" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plot_expansion_dataframe(total_expansion, \"Expansion capacity [GW]\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Plot of Installed capacity" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plot_expansion_dataframe(total_installed, \"Installed capacity [GW]\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Plot of optimal capacity" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plot_expansion_dataframe(total_optimal, \"Optimal capacity [GW]\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Plot of energy balance" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plot_expansion_dataframe(total_balance, \"Energy Balance in TWh\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "total_balance" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Check installed capacity in following year corresponds to optimal in the previous step" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "reorder_vector = [int(i/2) if i % 2 == 0 else len(total_optimal.columns) + int((i-1)/2) for i in range(2*len(total_optimal.columns))]\n", + "column_order = pd.concat([total_installed.columns.to_series(), total_optimal.columns.to_series()]).iloc[reorder_vector]\n", + "\n", + "pd.concat([total_optimal, total_installed], axis=1).loc[:, column_order]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Get retiring assets by planning_horizon" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# read costs from pypsa-earth\n", + "costs = pd.read_csv(root_myopic + f\"pypsa-earth/resources/costs.csv\", index_col=\"technology\")\n", + "\n", + "lifetime = costs.query(\"parameter == 'lifetime'\").value.rename(\"lifetime\")\n", + "lifetime[\"gas\"] = lifetime[\"OCGT\"] # gas tech is missing in costs.csv\n", + "\n", + "total_retiring_list = []\n", + "\n", + "# loop over the planning horizons to calculate the retiring capacity; skip the first planning horizon\n", + "for (i_pre, i_current) in zip(planning_horizons[:-1], planning_horizons[1:]):\n", + " # read previous network\n", + " n_pre = pypsa.Network(get_myopic_network(i_pre))\n", + "\n", + " subset_gen = n_pre.generators[n_pre.generators.carrier.isin(lifetime.index)]\n", + "\n", + " p_retiring = (\n", + " subset_gen.groupby(\"carrier\")\n", + " .apply(\n", + " lambda x: x.query(\n", + " f\"build_year < {i_current - lifetime[x.name.split('-')[0]]}\"\n", + " ).p_nom.sum()\n", + " )\n", + " .rename(f\"retiring - {i_current}\")\n", + " ).div(1e3).mul(-1.) # change sign to show negative values\n", + "\n", + " total_retiring_list.append(p_retiring)\n", + "\n", + "total_retiring = pd.concat(total_retiring_list, axis=1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plot_expansion_dataframe(total_retiring, \"Retired capacity [GW]\")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "pypsa-earth", + "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.13" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}