DR_potential_clustering.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Technical demand response potential clustering\n",
    "\n",
    "**Puporse and background**:\n",
    "* This notebook serves for clustering the technical demand response potentials which were collected in a previous meta-analysis (Kochems 2020).\n",
    "* A clustering routine is necessary since depicting all different units in the electricity market model would be computationally too expensive.\n",
    "\n",
    "## Method applied\n",
    "In the following, a brief description of the method applied is given for the sake of replicability.\n",
    "\n",
    "### Filtering of demand response categories\n",
    "The term \"demand response categories\" is introduced to describe the heterogeneous potential segmentation \n",
    "routines used in the publications evaluated within the meta-analysis. These demand response categories include \n",
    "processes, applications as well as entrire branches and mixtures of these categories.\n",
    "* In the first step, data for entire branches is filtered out since in most cases there are too few data points and further information on availability is lacking. Thus, only data for processes and applications remains.\n",
    "* In addition to that, categories covering different branches / appliances are filtered out since these would cause redundancies to publications analyzing the appliances in detail.\n",
    "\n",
    "### Prepare and manipulate data for further usage (clustering / modeling)\n",
    "In the second step, the data is manipulated for further usage in the clustering process as well as for the ensuing power market model analyses. This in turn consists of a few procedures:\n",
    "* The data is combined to an overall data set and missing values are interpolated using other parameter values as proxy, the median value per sector or 0s to fill nan values. The parameters are filtered in order only to include those which are needed in demand response modeling in a power market model.\n",
    "* A pairwhise correlation analysis using pearson's correlation coefficient is carried out in order to identify which parameters can be expressed through other ones since they show a high correlation. \n",
    "* Data is interpolated in order to remove inplausibilities. The data for the status quo is kept. The values for 2030 and 2050 are used to define trends in potential development. A linear interpolation is made in between. Also relations between positive and negative shifting potential as well as between shedding potential and installed capacity are kept. For the variable cost values, mean values are assigned. For investment expenses, the values for the status quo as well as the minimum value or at least 10% of the original value are kept. The minimum value is assigned to 2050 and a linear interpolation is made in between assuming a continued decline in costs. Fixed costs are assumed to be 2% of investment expenses.\n",
    "\n",
    "### Clustering of demand response categories\n",
    "A clustering of demand response categories is carried out in the second step. \n",
    "A k-means clustering approach is used (as an alternative, it is possible to choose agglomerative clustering using ward linkage).<br>\n",
    "Demand response categories are clustered using the (median values of the) following parameters (see also Steurer 2017, p. 83):\n",
    "* shifting duration\n",
    "* positive interference duration (shedding duration),\n",
    "* variable costs,\n",
    "* fixed costs and\n",
    "* specific investments.\n",
    "\n",
    "Some further aspects are worth mentioning:\n",
    "* Negative interference duration is not taken into account because some processes are only eligible for load shedding and hence don't have a negative interference duration. In addition to that, a strong correlation between positive and negative interference duration has been detected.\n",
    "* The clustering does not need to take into account the lower, middle and upper value for each parameter. A strong correlation between the values was determined which is why only the median values are used for the clustering. \n",
    "* Furthermore, the clustering is only carried out for the status quo and does not take development projections into account.\n",
    "* The distinction between different sectors an between shifting and shedding eligibilities is kept. Some heating and cooling applications for tcs and households are combined since they comprise basically the same technology and creating identical clusters would not make much sense.\n",
    "* For the aggregation of demand response parameters after clustering, a weighting by the available shifting resp. shedding capacity is carried out.\n",
    "\n",
    "### Determination of availability\n",
    "Since demand response potentials are time-dependent, availability has to be taken into account.<br>\n",
    "For the analysis, the individual availability time series in hourly resolution of the original \n",
    "demand response categories are put together by calculating capacity weighted averages for the identified \n",
    "demand response clusters.\n",
    "\n",
    "The availability time series are put together based on literature assumptions:\n",
    "* The largest amount of the availability time series for individual demand response categories were created within three bachelor theses based on literature assumptions (Benz 2019, Odeh 2019, Stange 2019). They were put together in a separate preprocessing jupyter notebook. The data output of this notebook in turn is read in here to form availability time series of demand response clusters.\n",
    "* Some processes haven't been analyzed in the bachelor theses resp. the literature. For these, either existing availability time series of very similar demand response categories are assigned or own assumptions are made."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Package imports\n",
    "* Standard imports: scikit-learn (sklearn) is used for the clustering since it has built-in clustering routines, such as K-means\n",
    "* User-defined functions:\n",
    "    * *group_potential*: Does a grouping of the clusters determined using given aggregation functions per parameter\n",
    "    * *write_multiple_sheets*: Used to write multiple DataFrames as sheets at once into an Excel workbook\n",
    "    * *map_column_names*: Maps column names of availability time series to the potential data column names.\n",
    "    * *determine_missing_columns*: Lists the columns for which availability time series information is lacking and assumptions are needed.\n",
    "    * *create_synthetic_profile_factors*: Creates synthetic profiles based on hourly, daily and monthly patterns\n",
    "    * *assign_availability_remaining*: Assigns availability time series for the remaining categories\n",
    "    * *get_top_abs_correlations*: Determines the strongest correlation within a given correlation matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": false
    }
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from sklearn.cluster import KMeans, AgglomerativeClustering\n",
    "from IPython.display import display\n",
    "import warnings\n",
    "\n",
    "from drpotentials.clustering_funcs import (\n",
    "    group_potential,\n",
    "    write_multiple_sheets,\n",
    "    map_column_names,\n",
    "    determine_missing_cols,\n",
    "    create_synthetic_profile_factors,\n",
    "    assign_availability_remaining,\n",
    "    get_top_abs_correlations\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": false
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "pd.options.mode.chained_assignment = None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Parameter settings\n",
    "Set path folder and filenames for reading in data.\n",
    "\n",
    "Further parameters for controlling the workflow:\n",
    "* *write_outputs*: If True, outputs, i.e. demand response parameterizations resulting from the meta-analysis will be written into Excel workbooks\n",
    "* *write_categories*: If True, the remaining categories will be written to Excel in order to match them with the\n",
    "availability data.\n",
    "* *adjust_potentials*: If True, availability time series information will be used to adjust the potential information. (I. e. if at max. 0.8 is reached, max potential will be set to 0.8 * max_potential)\n",
    "\n",
    "Parameters for controlling the clustering routines:\n",
    "* *quantile_cols*: Determine, which quantiles of the statistics DataFrame shall be used for demand response parameterization; Usually, 5%, 95% quantile as well as median are used. _**CAUTION**: If changes are applied, the order 'pessimistic', 'neutral', 'optimistic' must be kept, since it is used in the latter course of the analysis._\n",
    "* *cluster_parameters*: Determine, which demand response parameters to use for the clustering process. By default, these are the shifting duration, positive interference duration, variable costs as well as fixed costs.\n",
    "* *cluster_algo*: The clustering algorithmn to be used (\"KMeans\" or \"ward\")\n",
    "* *share_clusters*: Decide, how strong the original data will be reduced by giving a percentage of the original length. The cluster number is determined by the next higher integer. (Only applicable for k-means)\n",
    "* *distance_threshold*: Decide, what distance threshold shall be set for the hierarchical clustering using ward linkage, i.e., when the algorithm should terminate.\n",
    "* *print_clusters*: If True, prints out the clusters created (DataFrames)\n",
    "* *use_ava_ts_for_profiles*: If True, availability time series in positive direction will be directly used to derive load profiles for the demand response categories resp. clusters, else profiles from the demand regio disaggregator will be applied\n",
    "* *inflation_rate*: inflation rate is used to convert real costs values (as of 2020) to nominal terms which are used in the model _POMMES_."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": false
    }
   },
   "outputs": [],
   "source": [
    "# Set path folder(s) and filename(s) for reading in / writing data\n",
    "path_folder_in = \"./inputs/\"\n",
    "path_folder_stats = \"./out/stats/\"\n",
    "path_folder_availability = \"./out/availability/\"\n",
    "path_folder_parameterization = \"./out/parameterization/\"\n",
    "path_folder_plots = \"./out/plots/\"\n",
    "\n",
    "filename_in = \"Potenziale_Lastmanagement.xlsx\"\n",
    "filename_eligibility_in = \"eligibility_stats.csv\"\n",
    "filename_availability_in = \"availability_timeseries.xlsx\"\n",
    "filename_out = \"parameterization\"\n",
    "filename_corr_out = \"correlation\"\n",
    "filename_availability_out = \"availability_timeseries_clusters.csv\"\n",
    "filename_load_profiles_out = \"load_profile_timeseries_clusters\"\n",
    "\n",
    "# Set further parameters for controlling the workflow\n",
    "write_outputs = True\n",
    "write_categories = True\n",
    "adjust_potentials = True\n",
    "\n",
    "# Determine clustering approach\n",
    "quantile_cols = [\"5%\", \"50%\", \"95%\"]\n",
    "cluster_parameters = [\n",
    "    \"shifting_duration\", \"interference_duration_pos\", \n",
    "    \"variable_costs\", \"fixed_costs\", \"specific_investments\"\n",
    "]\n",
    "cluster_parameters = [a + \"_\" + b for a in quantile_cols for b in cluster_parameters]\n",
    "\n",
    "cluster_algo = \"KMeans\"\n",
    "share_clusters = 0.1\n",
    "distance_threshold = 1000\n",
    "print_clusters = True\n",
    "use_ava_ts_for_profiles = False\n",
    "inflation_rate = 0.02"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SMALL_SIZE = 12\n",
    "MEDIUM_SIZE = 14\n",
    "BIGGER_SIZE = 15\n",
    "\n",
    "plt.rc('font', size=SMALL_SIZE)          # controls default text sizes\n",
    "plt.rc('axes', titlesize=BIGGER_SIZE)    # fontsize of the axes title\n",
    "plt.rc('axes', labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels\n",
    "plt.rc('xtick', labelsize=SMALL_SIZE)    # fontsize of the tick labels\n",
    "plt.rc('ytick', labelsize=SMALL_SIZE)    # fontsize of the tick labels\n",
    "plt.rc('legend', fontsize=SMALL_SIZE)    # legend fontsize\n",
    "plt.rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure title"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Read in and filter data\n",
    "* Read in the categories data and filter out branches as well as power-to-X-technologies which won't be considered anymore.\n",
    "* Read in the stats information on the demand response parameters from the previous meta-analysis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in and filter demand response categories\n",
    "* Read in demand response categories\n",
    "* Drop entire branches as well as categories conflicting with others or outside of scope (Power-to-X other than Power-to-Heat for space heating).\n",
    "* Show the original number of categories and print the remaining number after filtering (without duplicates since categories may be used within different sectors)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": false
    }
   },
   "outputs": [],
   "source": [
    "categories = pd.read_excel(path_folder_in+filename_in, sheet_name=\"Kategorien_neu\", index_col=0)\n",
    "colors = pd.read_excel(path_folder_in+filename_in, sheet_name=\"Farben\", index_col=0)\n",
    "\n",
    "print(f\"Number of original categories:\\t{categories.shape[0]}\")\n",
    "\n",
    "categories = categories[\n",
    "    categories[\"Nutzung?\"] == 1 \n",
    "    & ~categories[\"Einstufung\"].isin([\"Branche\", \"Power-to-X\"])\n",
    "]\n",
    "\n",
    "categories = categories.drop_duplicates(subset=\"Prozesskategorie\")\n",
    "categories = categories.set_index([\"Prozesskategorie\"], drop=True)\n",
    "\n",
    "print(f\"Number of remaining categories:\\t{categories.shape[0]}\")\n",
    "\n",
    "# Show the remaining demand response categories which are evaluated\n",
    "# list(categories.index.values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in demand response parameters data\n",
    "* Assign each demand response parameter an aggregation function to be used after clustering (sum or mean).\n",
    "* Determine for which parameters to swap the order of preference.\n",
    "> _Note: This is necessary, because in some cases minimum values are needed for an optimistic \n",
    "demand response projection and maximum for a pessimistic one, e.g. for minimum load factor. <br>\n",
    "Hence, for these parameters, min is exchanged for max etc._\n",
    "* Determine sectors and years for which parameter information existis.\n",
    "* Determine, which categories to join for household and tcs sector due to the same appliances being used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": false
    }
   },
   "outputs": [],
   "source": [
    "# Assgin each parameter the aggregation function to be used\n",
    "parameters_agg_dict = {\n",
    "    \"activation_duration\": \"mean\", \n",
    "    \"ave_load\": \"mean\", \n",
    "    \"fixed_costs\": \"mean\", \n",
    "    \"installed_cap\": \"sum\",\n",
    "    \"interference_duration_neg\": \"mean\", \n",
    "    \"interference_duration_pos\": \"mean\",\n",
    "    \"interference_duration_pos_shed\": \"mean\",\n",
    "    \"max_load\": \"mean\", \n",
    "    \"maximum_activations_year\": \"mean\", \n",
    "    \"maximum_activations_year_shed\": \"mean\",\n",
    "    \"min_load\": \"mean\", \n",
    "    \"potential_neg_overall\": \"sum\",\n",
    "    \"potential_pos_overall\": \"sum\",\n",
    "    \"potential_pos_overall_shed\": \"sum\",\n",
    "    \"regeneration_duration\": \"mean\", \n",
    "    \"shiftable_share\": \"mean\", \n",
    "    \"shifting_duration\": \"mean\",\n",
    "    \"specific_investments\": \"mean\", \n",
    "    \"variable_costs\": \"mean\",\n",
    "    \"variable_costs_shed\": \"mean\"\n",
    "}\n",
    "\n",
    "# Determine for each parameter whether or not to swap values\n",
    "parameters_swap_dict = {\n",
    "    \"activation_duration\": True, \n",
    "    \"ave_load\": True, \n",
    "    \"fixed_costs\": True, \n",
    "    \"installed_cap\": False,\n",
    "    \"interference_duration_neg\": False, \n",
    "    \"interference_duration_pos\": False,\n",
    "    \"interference_duration_pos_shed\": False,\n",
    "    \"max_load\": False, \n",
    "    \"maximum_activations_year\": False, \n",
    "    \"maximum_activations_year_shed\": False,\n",
    "    \"min_load\": True, \n",
    "    \"potential_neg_overall\": False,\n",
    "    \"potential_pos_overall\": False,\n",
    "    \"potential_pos_overall_shed\": False,\n",
    "    \"regeneration_duration\": True, \n",
    "    \"shiftable_share\": False, \n",
    "    \"shifting_duration\": False,\n",
    "    \"specific_investments\": True, \n",
    "    \"variable_costs\": True,\n",
    "    \"variable_costs_shed\": True  \n",
    "}\n",
    "\n",
    "# Map columns for swapping\n",
    "swap_cols = {\n",
    "    \"min\": \"max\",\n",
    "    \"5%\": \"95%\",\n",
    "    \"10%\": \"90%\",\n",
    "    \"25%\": \"75%\"\n",
    "}\n",
    "\n",
    "sectors = [\"ind\", \"tcs\", \"hoho\"]\n",
    "\n",
    "years = [\"SQ\", \"2020\", \"2025\", \"2030\", \"2035\", \"2040\", \"2045\", \"2050\"]\n",
    "\n",
    "to_join = [\"Nachtspeicherheizungen\", \"Warmwasserbereitstellung\", \"Wärmepumpen\", \"Klimakälte\"]\n",
    "to_drop = {\"Prozesskälte\": \"hoho\"}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a common data basis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Count the number of params for which no data is available\n",
    "count_ignored = 0\n",
    "\n",
    "to_concat = []\n",
    "for parameter, swap_param in parameters_swap_dict.items():\n",
    "    for year in years:\n",
    "        for sector in sectors:\n",
    "            try:\n",
    "                new_df = pd.read_csv(\n",
    "                    path_folder_stats+parameter+\"_\"+sector+\"_stats_\" + year + \".csv\",\n",
    "                    sep=\";\",\n",
    "                    index_col=0\n",
    "                ).T\n",
    "                new_df.index.name = \"Prozesskategorie\"\n",
    "\n",
    "                # Change the order of appearance if swap is needed for parameter\n",
    "                if swap_param:\n",
    "                    for k, v in swap_cols.items():\n",
    "                        new_df[k + \"_copy\"] = new_df[k]\n",
    "                        new_df[v + \"_copy\"] = new_df[v]\n",
    "                        new_df[k] = new_df[v + \"_copy\"]\n",
    "                        new_df[v] = new_df[k + \"_copy\"]\n",
    "                        new_df.drop(columns=[k + \"_copy\", v + \"_copy\"], inplace=True)\n",
    "                    \n",
    "                new_df[\"parameter\"] = parameter\n",
    "                new_df[\"sector\"] = sector\n",
    "                new_df[\"year\"] = year\n",
    "                to_concat.append(new_df)\n",
    "            except FileNotFoundError:\n",
    "                count_ignored += 1\n",
    "                continue\n",
    "\n",
    "# Put everything into one common pandas.DataFrame\n",
    "all_params_df = pd.concat(to_concat, sort=False)\n",
    "                \n",
    "print(f\"Overall number of params (sectors, years, params): \"\n",
    "      f\"{len(parameters_agg_dict) * (len(sectors)) * len(years)}\")\n",
    "print(f\"Number of params not eligible for evaluation: {count_ignored}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Drop values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Filter out the categories to be used and drop certain ones for dedicated sector only\n",
    "all_params_df = all_params_df.loc[all_params_df.index.isin(categories.index)]\n",
    "all_params_df.set_index(\"sector\", append=True, inplace=True)\n",
    "all_params_df.drop(index={(k, v) for k, v in to_drop.items()}, inplace=True)\n",
    "all_params_df.reset_index(level=1, inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Isolate an uncombined version of the overall data set\n",
    "\n",
    "> _Note:_\n",
    "> * This does not combine space climate and heating processes for tcs and households.\n",
    "> * The uncombined version is needed because differing profiles for tcs / hoho have to be assigned in the later course."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_params_df_unjoined = all_params_df.copy()\n",
    "all_params_df_unjoined = all_params_df.set_index([\"year\", \"parameter\", \"sector\"], append=True).round(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Combine duplicates for tcs & hoho\n",
    "> _Note:_\n",
    "> * This combined version is needed since space climate and heating processes for tcs and households shall not be treated differently in the clustering routine that follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extend index and filter for duplicate values\n",
    "all_params_df.set_index([\"year\", \"parameter\"], append=True, inplace=True)\n",
    "filter_duplicates = (\n",
    "    all_params_df.index.duplicated(keep=False)\n",
    "    & all_params_df.sector.isin([\"hoho\", \"tcs\"])\n",
    ")\n",
    "\n",
    "# Get duplicates within duplicates (i.e. filter out the values where \n",
    "# the same demand response categories occur in tcs and industry sector)\n",
    "duplicates_df = all_params_df.loc[filter_duplicates]\n",
    "duplicates_df = duplicates_df[duplicates_df.index.duplicated(keep=False)]\n",
    "\n",
    "# Filter out the remaining duplicate values from original DataFrame\n",
    "keys = list(duplicates_df.columns.values)\n",
    "i1 = all_params_df.set_index(keys).index\n",
    "i2 = duplicates_df.set_index(keys).index\n",
    "no_duplicates_df = all_params_df.loc[~i1.isin(i2)]\n",
    "\n",
    "# Assign certain demand response categories to a combined tcs & hoho sector (\"tcs+hoho\")\n",
    "for el in to_join:\n",
    "    if el in no_duplicates_df.index:\n",
    "        no_duplicates_df.loc[:,\"help_sector\"] = np.where(\n",
    "            no_duplicates_df.loc[:,\"sector\"].values == \"ind\", \"ind\", \"tcs+hoho\"\n",
    "        )\n",
    "        no_duplicates_df.loc[el, \"sector\"] = no_duplicates_df.loc[el, \"help_sector\"].values\n",
    "        no_duplicates_df = no_duplicates_df.drop([\"help_sector\"], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Drop index for grouping\n",
    "duplicates_df.reset_index(drop=False, inplace=True)\n",
    "\n",
    "# Introduce DataFrame for grouping duplicates (sector \"tcs+hoho\")\n",
    "grouped_duplicates_df = pd.DataFrame()\n",
    "\n",
    "for parameter, param_agg_rule in parameters_agg_dict.items():\n",
    "    param_duplicates_df = duplicates_df.loc[duplicates_df[\"parameter\"] == parameter]\n",
    "    if not param_duplicates_df.empty:\n",
    "        param_duplicates_df = param_duplicates_df.groupby(\n",
    "            [\"Prozesskategorie\", \"year\", \"parameter\"]).agg({\n",
    "                \"count\": \"sum\",\n",
    "                \"min\": param_agg_rule,\n",
    "                \"5%\": param_agg_rule,\n",
    "                \"10%\": param_agg_rule,\n",
    "                \"25%\": param_agg_rule,\n",
    "                \"50%\": param_agg_rule,\n",
    "                \"mean\": param_agg_rule,\n",
    "                \"std\": param_agg_rule,\n",
    "                \"75%\": param_agg_rule,\n",
    "                \"90%\": param_agg_rule,\n",
    "                \"95%\": param_agg_rule,\n",
    "                \"max\": param_agg_rule,\n",
    "            })\n",
    "        grouped_duplicates_df = pd.concat([grouped_duplicates_df, param_duplicates_df])\n",
    "        \n",
    "grouped_duplicates_df[\"sector\"] = \"tcs+hoho\"\n",
    "\n",
    "# Combine the non-duplicated data with the consolidated tcs+hoho data\n",
    "all_params_df = pd.concat([no_duplicates_df, grouped_duplicates_df], sort=False)\n",
    "\n",
    "# Create unique index and round to four digits\n",
    "all_params_df.set_index(\"sector\", append=True, inplace=True)\n",
    "all_params_df = all_params_df.round(4)\n",
    "\n",
    "sectors.append(\"tcs+hoho\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in and append eligibility information\n",
    "* Read in eligibility information extracted from potential evaluation notebook\n",
    "* Assign eligibility for shifting / shedding when more than 50% of the data sources indicate elibility to do so\n",
    "* Combine eligibility with demand response parameters data set and filter out processes found to be not eligible"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eligibility = pd.read_csv(path_folder_in+filename_eligibility_in, index_col=0)\n",
    "eligibility = eligibility[eligibility.index.isin(categories.index)]\n",
    "\n",
    "# Determine relative values and assign binaries for shifting / shedding eligibility\n",
    "eligibility = (\n",
    "    eligibility[[\"Lastverschiebung\", \"Lastverzicht\"]].div(eligibility[\"Anzahl\"], axis=0)\n",
    ")\n",
    "eligibility[\"shifting\"] = np.where(eligibility[\"Lastverschiebung\"] >= 0.5, 1, 0)\n",
    "eligibility[\"shedding\"] = np.where(eligibility[\"Lastverzicht\"] >= 0.5, 1, 0)\n",
    "\n",
    "# Introduce strings to distinct different shifting / shedding eligibilities\n",
    "conditions = [\n",
    "    (eligibility[\"shifting\"] == 1) & (eligibility[\"shedding\"] == 1),\n",
    "    (eligibility[\"shifting\"] == 1) & (eligibility[\"shedding\"] != 1),\n",
    "    (eligibility[\"shifting\"] != 1) & (eligibility[\"shedding\"] == 1)\n",
    "]\n",
    "choices = [\n",
    "    \"shift_shed\", \"shift_only\", \"shed_only\"\n",
    "]\n",
    "eligibility[\"kind\"] = np.select(\n",
    "    conditions, choices, \"not_eligible\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Combine eligibility data with parameter data set - combined data set\n",
    "all_params_df = all_params_df.reset_index(level=[1, 2, 3])\n",
    "all_params_df[\"kind\"] = eligibility[\"kind\"]\n",
    "all_params_df = all_params_df.loc[all_params_df[\"kind\"] != \"not_eligible\"].set_index(\n",
    "    [\"year\", \"parameter\", \"sector\", \"kind\"], append=True\n",
    ")\n",
    "\n",
    "# Uncombined data set\n",
    "all_params_df_unjoined = all_params_df_unjoined.reset_index(level=[1, 2, 3])\n",
    "all_params_df_unjoined[\"kind\"] = eligibility[\"kind\"]\n",
    "all_params_df_unjoined = all_params_df_unjoined.loc[all_params_df_unjoined[\"kind\"] != \"not_eligible\"].set_index(\n",
    "    [\"year\", \"parameter\", \"sector\", \"kind\"], append=True\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Prepare and manipulate data for further usage"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reshape, fill nan values and re-combine data for cluster analysis\n",
    "* Reshape data such that parameters are columns and process categories form the index.\n",
    "* Fillna values based on defined rules:\n",
    "    * If positive or negative potential is set while the other is not, define negative potential to equal positive one (or vice versa; assumption of symmetrical potentials).\n",
    "    * Use positive or negative potential as a proxy where installed capacity is missing.\n",
    "    * Assign median values per sector and year for parameters dependent on an attribution to positive, negative or shedding potential\n",
    "    * Assign 0 values to all remaining missing values.\n",
    "    * Fill stll remaining Nan values with zero entries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# All parameter columns\n",
    "all_param_cols = [\n",
    "    (a, b) for a in quantile_cols for b in parameters_agg_dict.keys()\n",
    "]\n",
    "\n",
    "# Define potential columns\n",
    "potential_neg_cols = [\n",
    "    (a, b) for a in quantile_cols for b in [\"potential_neg_overall\"]\n",
    "]\n",
    "potential_pos_cols = [\n",
    "    (a, b) for a in quantile_cols for b in [\"potential_pos_overall\"]\n",
    "]\n",
    "potential_pos_shed_cols = [\n",
    "    (a, b) for a in quantile_cols for b in [\"potential_pos_overall_shed\"]\n",
    "]\n",
    "installed_cap_cols = [\n",
    "    (a, b) for a in quantile_cols for b in [\"installed_cap\"]\n",
    "]\n",
    "\n",
    "# Fill nan for parameters when any of potential columns is not None\n",
    "cols_any = [\n",
    "    \"activation_duration\",\n",
    "    \"ave_load\",\n",
    "    \"fixed_costs\",\n",
    "    \"installed_cap\",\n",
    "    \"max_load\",\n",
    "    \"min_load\",\n",
    "    \"regeneration_duration\",\n",
    "    \"specific_investments\"\n",
    "]\n",
    "cols_any = [(a, b) for a in quantile_cols for b in cols_any]\n",
    "\n",
    "# Fill nan for parameters when shifting is possible\n",
    "cols_shift = [\n",
    "    \"maximum_activations_year\",\n",
    "    \"shiftable_share\",\n",
    "    \"shifting_duration\",\n",
    "    \"variable_costs\"\n",
    "]\n",
    "cols_shift = [(a, b) for a in quantile_cols for b in cols_shift]\n",
    "\n",
    "shiftable_share_cols = [\n",
    "    (a, b) for a in quantile_cols for b in [\"shiftable_share\"]\n",
    "]\n",
    "\n",
    "# Fill nan for parameters when negative / positive / shedding potential information is given\n",
    "cols_neg = [\n",
    "    \"interference_duration_neg\",\n",
    "]\n",
    "cols_neg = [(a, b) for a in quantile_cols for b in cols_neg]\n",
    "\n",
    "cols_pos = [\n",
    "    \"interference_duration_pos\"\n",
    "]\n",
    "cols_pos = [(a, b) for a in quantile_cols for b in cols_pos]\n",
    "\n",
    "cols_shed = [\n",
    "    \"interference_duration_pos_shed\",\n",
    "    \"maximum_activations_year_shed\",\n",
    "    \"variable_costs_shed\"\n",
    "]\n",
    "cols_shed = [(a, b) for a in quantile_cols for b in cols_shed]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_params_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reshape data for correlation analysis and clustering - combined data set\n",
    "all_params_reshaped = pd.DataFrame()\n",
    "\n",
    "for year in years:\n",
    "\n",
    "    # Due to multiple indices (duplicate categories for ind / tcs) needs to be done by sector\n",
    "    for sector in sectors:\n",
    "        reshaped_df = all_params_df.loc[\n",
    "            (all_params_df.index.get_level_values(1) == year)\n",
    "            & (all_params_df.index.get_level_values(3) == sector)\n",
    "        ]\n",
    "\n",
    "        # Restructure the data: process categories row-wise, parameters column-wise\n",
    "        reshaped_df.drop(\n",
    "            columns=[\n",
    "                col for col in reshaped_df.columns \n",
    "                if col not in quantile_cols\n",
    "            ], \n",
    "            inplace=True\n",
    "        )\n",
    "\n",
    "        reshaped_df.reset_index(drop=False, inplace=True)\n",
    "        reshaped_df = reshaped_df.pivot(\n",
    "            index=[\"Prozesskategorie\", \"sector\", \"year\", \"kind\"], \n",
    "            columns=[\"parameter\"]\n",
    "        )\n",
    "        reshaped_df = reshaped_df.reindex(columns=all_param_cols)\n",
    "\n",
    "        \n",
    "        # Determine processes, for which (no) potential / installed capacity information is given\n",
    "        neg_nan_idx = reshaped_df[\n",
    "            reshaped_df[[col for col in potential_neg_cols]].isna().all(axis=1)\n",
    "        ].index\n",
    "\n",
    "        neg_not_nan_idx = [\n",
    "            idx for idx in reshaped_df.index if idx not in neg_nan_idx \n",
    "        ]\n",
    "\n",
    "        pos_nan_idx = reshaped_df[\n",
    "            reshaped_df[[col for col in potential_pos_cols]].isna().all(axis=1)\n",
    "        ].index\n",
    "\n",
    "        pos_not_nan_idx = [\n",
    "            idx for idx in reshaped_df.index if idx not in pos_nan_idx\n",
    "        ]\n",
    "\n",
    "        pos_shed_nan_idx = reshaped_df[\n",
    "            reshaped_df[[col for col in potential_pos_shed_cols]].isna().all(axis=1)\n",
    "        ].index\n",
    "\n",
    "        pos_shed_not_nan_idx = [\n",
    "            idx for idx in reshaped_df.index if idx not in pos_shed_nan_idx\n",
    "        ]\n",
    "        \n",
    "        installed_nan_idx = reshaped_df[\n",
    "            reshaped_df[[col for col in installed_cap_cols]].isna().all(axis=1)\n",
    "        ].index\n",
    "\n",
    "        neg_missing = neg_nan_idx.difference(pos_nan_idx)\n",
    "        pos_missing = pos_nan_idx.difference(neg_nan_idx)\n",
    "        installed_cap_missing_fill_pos = installed_nan_idx.difference(pos_nan_idx)\n",
    "        installed_cap_missing_fill_neg = installed_nan_idx.difference(neg_nan_idx)\n",
    "        \n",
    "        # Fill categories missing negative potential with positive one as a proxy\n",
    "        reshaped_df.loc[neg_missing, potential_neg_cols] = (\n",
    "            reshaped_df.loc[neg_missing, potential_pos_cols].values\n",
    "        )\n",
    "\n",
    "        # Fill categories missing positive potential with negative one as a proxy\n",
    "        reshaped_df.loc[pos_missing, potential_pos_cols] = (\n",
    "            reshaped_df.loc[pos_missing, potential_neg_cols].values\n",
    "        )\n",
    "        \n",
    "        # Fill categories missing installed capacity with upper estimate \n",
    "        # for positive or negative potential as a proxy\n",
    "        try:\n",
    "            reshaped_df.loc[installed_cap_missing_fill_pos, installed_cap_cols] = (\n",
    "                reshaped_df.loc[installed_cap_missing_fill_pos, potential_pos_cols[-1]]\n",
    "            )\n",
    "        except ValueError:\n",
    "            pass\n",
    "        try:\n",
    "            reshaped_df.loc[installed_cap_missing_fill_neg, installed_cap_cols] = (\n",
    "                reshaped_df.loc[installed_cap_missing_fill_neg, potential_pos_cols[-1]]\n",
    "            )\n",
    "        except ValueError:\n",
    "            pass\n",
    "\n",
    "        # Fill non-potential parameters dependent on potentials given\n",
    "        for col in cols_neg + cols_shift + cols_any:\n",
    "            reshaped_df.loc[neg_not_nan_idx, col] = (\n",
    "                reshaped_df.loc[neg_not_nan_idx, col].fillna(\n",
    "                    reshaped_df[col].median()\n",
    "                )\n",
    "            )\n",
    "\n",
    "        for col in cols_pos + cols_shift + cols_any:\n",
    "            reshaped_df.loc[pos_not_nan_idx, col] = (\n",
    "                reshaped_df.loc[pos_not_nan_idx, col].fillna(\n",
    "                    reshaped_df[col].median() \n",
    "                )\n",
    "            )\n",
    "\n",
    "        for col in cols_shed + cols_any:\n",
    "            reshaped_df.loc[pos_shed_not_nan_idx, col] = (\n",
    "                reshaped_df.loc[pos_shed_not_nan_idx, col].fillna(\n",
    "                    reshaped_df[col].median() \n",
    "                )\n",
    "            )\n",
    "        \n",
    "        # Combine\n",
    "        all_params_reshaped = pd.concat([all_params_reshaped, reshaped_df])\n",
    "    \n",
    "    # Drop processes that have neither potential nor installed capacity information\n",
    "    all_params_reshaped.drop(\n",
    "        index=all_params_reshaped.loc[\n",
    "            all_params_reshaped.loc[\n",
    "                :, \n",
    "                installed_cap_cols \n",
    "                + potential_pos_cols \n",
    "                + potential_neg_cols \n",
    "                + potential_pos_shed_cols\n",
    "            ].isna().all(axis=1)\n",
    "        ].index,\n",
    "        inplace=True\n",
    "    )\n",
    "    \n",
    "    # Fill remaining nans\n",
    "    all_params_reshaped[shiftable_share_cols] = all_params_reshaped[shiftable_share_cols].fillna(1)\n",
    "    all_params_reshaped = all_params_reshaped.fillna(0).round(2)\n",
    "    all_params_reshaped.sort_index(axis=1, level=1, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reshape data for clustering - uncombined data set\n",
    "all_params_reshaped_unjoined = pd.DataFrame()\n",
    "\n",
    "for year in years:\n",
    "\n",
    "    # Due to multiple indices (duplicate categories for ind / tcs) needs to be done by sector\n",
    "    for sector in sectors:\n",
    "        reshaped_df = all_params_df_unjoined.loc[\n",
    "            (all_params_df_unjoined.index.get_level_values(1) == year)\n",
    "            & (all_params_df_unjoined.index.get_level_values(3) == sector)\n",
    "        ]\n",
    "\n",
    "        # Restructure the data: process categories row-wise, parameters column-wise\n",
    "        all_params_df_unjoined.drop(\n",
    "            columns=[\n",
    "                col for col in reshaped_df.columns \n",
    "                if col not in quantile_cols\n",
    "            ], \n",
    "            inplace=True\n",
    "        )\n",
    "\n",
    "        reshaped_df.reset_index(drop=False, inplace=True)\n",
    "        reshaped_df = reshaped_df.pivot(\n",
    "            index=[\"Prozesskategorie\", \"sector\", \"year\", \"kind\"], \n",
    "            columns=[\"parameter\"]\n",
    "        )\n",
    "        reshaped_df = reshaped_df.reindex(columns=all_param_cols)\n",
    "        \n",
    "        # Determine processes, for which (no) potential information is given\n",
    "        neg_nan_idx = reshaped_df[\n",
    "            reshaped_df[[col for col in potential_neg_cols]].isna().all(axis=1)\n",
    "        ].index\n",
    "\n",
    "        neg_not_nan_idx = [\n",
    "            idx for idx in reshaped_df.index if idx not in neg_nan_idx \n",
    "        ]\n",
    "\n",
    "        pos_nan_idx = reshaped_df[\n",
    "            reshaped_df[[col for col in potential_pos_cols]].isna().all(axis=1)\n",
    "        ].index\n",
    "\n",
    "        pos_not_nan_idx = [\n",
    "            idx for idx in reshaped_df.index if idx not in pos_nan_idx\n",
    "        ]\n",
    "\n",
    "        pos_shed_nan_idx = reshaped_df[\n",
    "            reshaped_df[[col for col in potential_pos_shed_cols]].isna().all(axis=1)\n",
    "        ].index\n",
    "\n",
    "        pos_shed_not_nan_idx = [\n",
    "            idx for idx in reshaped_df.index if idx not in pos_shed_nan_idx\n",
    "        ]\n",
    "\n",
    "        installed_nan_idx = reshaped_df[\n",
    "            reshaped_df[[col for col in installed_cap_cols]].isna().all(axis=1)\n",
    "        ].index\n",
    "\n",
    "        neg_missing = neg_nan_idx.difference(pos_nan_idx)\n",
    "        pos_missing = pos_nan_idx.difference(neg_nan_idx)\n",
    "        installed_cap_missing_fill_pos = installed_nan_idx.difference(pos_nan_idx)\n",
    "        installed_cap_missing_fill_neg = installed_nan_idx.difference(neg_nan_idx)\n",
    "\n",
    "        # Fill categories missing negative potential with positive one as a proxy\n",
    "        reshaped_df.loc[neg_missing, potential_neg_cols] = (\n",
    "            reshaped_df.loc[neg_missing, potential_pos_cols].values\n",
    "        )\n",
    "\n",
    "        # Fill categories missing positive potential with negative one as a proxy\n",
    "        reshaped_df.loc[pos_missing, potential_pos_cols] = (\n",
    "            reshaped_df.loc[pos_missing, potential_neg_cols].values\n",
    "        )\n",
    "\n",
    "        # Fill categories missing installed capacity with upper estimate \n",
    "        # for positive or negative potential as a proxy\n",
    "        try:\n",
    "            reshaped_df.loc[installed_cap_missing_fill_pos, installed_cap_cols] = (\n",
    "                reshaped_df.loc[installed_cap_missing_fill_pos, potential_pos_cols[-1]]\n",
    "            )\n",
    "        except ValueError:\n",
    "            pass\n",
    "        try:\n",
    "            reshaped_df.loc[installed_cap_missing_fill_neg, installed_cap_cols] = (\n",
    "                reshaped_df.loc[installed_cap_missing_fill_neg, potential_pos_cols[-1]]\n",
    "            )\n",
    "        except ValueError:\n",
    "            pass\n",
    "\n",
    "        # Fill non-potential parameters dependent on potentials given\n",
    "        for col in cols_neg + cols_shift + cols_any:\n",
    "            reshaped_df.loc[neg_not_nan_idx, col] = (\n",
    "                reshaped_df.loc[neg_not_nan_idx, col].fillna(\n",
    "                    reshaped_df[col].median() \n",
    "                )\n",
    "            )\n",
    "\n",
    "        for col in cols_pos + cols_shift + cols_any:\n",
    "            reshaped_df.loc[pos_not_nan_idx, col] = (\n",
    "                reshaped_df.loc[pos_not_nan_idx, col].fillna(\n",
    "                    reshaped_df[col].median() \n",
    "                )\n",
    "            )\n",
    "\n",
    "        for col in cols_shed + cols_any:\n",
    "            reshaped_df.loc[pos_shed_not_nan_idx, col] = (\n",
    "                reshaped_df.loc[pos_shed_not_nan_idx, col].fillna(\n",
    "                    reshaped_df[col].median() \n",
    "                )\n",
    "            )\n",
    "        \n",
    "        # Combine\n",
    "        all_params_reshaped_unjoined = pd.concat([all_params_reshaped_unjoined, reshaped_df])\n",
    "    \n",
    "    # Drop processes that have neither potential nor installed capacity information\n",
    "    all_params_reshaped_unjoined.drop(\n",
    "        index=all_params_reshaped_unjoined.loc[\n",
    "            all_params_reshaped_unjoined.loc[\n",
    "                :, \n",
    "                installed_cap_cols \n",
    "                + potential_pos_cols \n",
    "                + potential_neg_cols \n",
    "                + potential_pos_shed_cols\n",
    "            ].isna().all(axis=1)\n",
    "        ].index,\n",
    "        inplace=True\n",
    "    )\n",
    "    \n",
    "    # Fill remaining nans\n",
    "    all_params_reshaped_unjoined[shiftable_share_cols] = all_params_reshaped_unjoined[shiftable_share_cols].fillna(1)\n",
    "    all_params_reshaped_unjoined = all_params_reshaped_unjoined.fillna(0).round(2)\n",
    "    all_params_reshaped_unjoined.sort_index(axis=1, level=1, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_params_reshaped.tail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_params_reshaped_unjoined.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Do a correlation analysis of parameters in order to derive cluster parameters\n",
    "General idea:\n",
    "* Check pairwise correlations between parameters to detect autocorrelation and derive parameters used for clustering.\n",
    "* If correlation between to parameters is above 0.8, choose one of both."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "corr_matrix = all_params_reshaped.loc[\n",
    "    all_params_reshaped.index.get_level_values(2) == \"SQ\"\n",
    "].corr().round(2).fillna(0)\n",
    "\n",
    "# Transform to 1-dimensional indices\n",
    "corr_matrix[\"new_index\"] = (\n",
    "    corr_matrix.index.get_level_values(0) + \"_\" + corr_matrix.index.get_level_values(1)\n",
    ")\n",
    "corr_matrix.set_index(\"new_index\", drop=True, inplace=True)\n",
    "corr_matrix.columns = (\n",
    "    corr_matrix.columns.get_level_values(0) + \"_\" + corr_matrix.columns.get_level_values(1)\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display the top 10 largest correlations\n",
    "print(\"Top Absolute Correlations\")\n",
    "print(75 * \"-\")\n",
    "print(get_top_abs_correlations(corr_matrix, n=10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Determine and store top correlations\n",
    "top_corr_series = get_top_abs_correlations(corr_matrix, threshold=0.8)\n",
    "top_corr_series.name = \"pearson_correlation\"\n",
    "corr_dict = {}\n",
    "\n",
    "# Filter for params and quantiles\n",
    "ix1_slice_end = top_corr_series.index.get_level_values(0).str.split(\n",
    "    \"_\", 1, expand=True\n",
    ").get_level_values(1)\n",
    "ix2_slice_end = top_corr_series.index.get_level_values(1).str.split(\n",
    "    \"_\", 1, expand=True\n",
    ").get_level_values(1)\n",
    "ix1_slice_start = top_corr_series.index.get_level_values(0).str.split(\n",
    "    \"_\", 1, expand=True\n",
    ").get_level_values(0)\n",
    "ix2_slice_start = top_corr_series.index.get_level_values(1).str.split(\n",
    "    \"_\", 1, expand=True\n",
    ").get_level_values(0)\n",
    "\n",
    "corr_dict[\"same_params\"] = top_corr_series[ix1_slice_end == ix2_slice_end]\n",
    "corr_dict[\"different_params\"] = top_corr_series[ix1_slice_end != ix2_slice_end]\n",
    "corr_dict[\"only_medians\"] = top_corr_series[\n",
    "    (ix1_slice_end != ix2_slice_end)\n",
    "    & (ix1_slice_start == \"50%\") \n",
    "    & (ix2_slice_start == \"50%\")\n",
    "]\n",
    "\n",
    "if write_outputs:\n",
    "    write_multiple_sheets(\n",
    "        corr_dict, path_folder_parameterization, filename_corr_out + \".xlsx\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Evaluate correlations for medians only\n",
    "ix_slice_start = corr_matrix.index.str.split(\n",
    "    \"_\", 1, expand=True\n",
    ").get_level_values(0)\n",
    "\n",
    "corr_matrix_50 = corr_matrix.loc[\n",
    "    ix_slice_start == \"50%\", \n",
    "    [col for col in corr_matrix if \"50%\" in col]\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show heatmap for all correlations\n",
    "fig, ax = plt.subplots(figsize=(20, 15))\n",
    "_ = sns.heatmap(corr_matrix, cmap=\"RdGy\", vmin=-1.0, vmax=1.0)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show heatmap for median values only\n",
    "fig, ax = plt.subplots(figsize=(20, 15))\n",
    "_ = sns.heatmap(corr_matrix_50, cmap=\"RdGy\", vmin=-1.0, vmax=1.0)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create consistent trends for future developments\n",
    "Overview on approach used:\n",
    "1. Introduce **linear interpolation** for capital cost-related potential parameters:\n",
    "    * Potential parameters: The values for the status quo, 2030 and 2050 are used. Linear interpolation is made in between.\n",
    "    * Cost parameters: For specific investments, the maximum of the minimum cost value and 10% of original costs is assigned to 2050. This is for two reasons:\n",
    "        * A cost degression over time is likely.\n",
    "        * Nonetheless, some studies containing extremely low values, such as 0 or close to that, assume that no additional investments is deemed necessary which is not in line with the research design applied here.\n",
    "2. For variable costs, calculate **averages** over all years and assume constant costs in real terms, thus inflationing values in nominal terms.\n",
    "3. **Assign remaining** parameters (mostly time-related ones) the same as for the status quo for every year."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Only a **subset of the parameters** needs to be further analyzed since not every parameter is needed for modelling:\n",
    "* Time-dependent parameters are assumed constant. These comprise:\n",
    "    * activation duration\n",
    "    * interference duration (both pos and neg) and shifting duration\n",
    "    * regeneration duration\n",
    "* Other parameters are not really resp. not directly used in the modeling approaches for DR. These comprise:\n",
    "    * average, minimum and maximum load\n",
    "* This leads to the following remaining parameters focussing on costs and potentials. Since the correlation analysis showed high redundandencies for the potential parameters, only the following remaining parameters will be further analyzed:\n",
    "    * potential positive overall\n",
    "    * potential positive overall for shedding\n",
    "    * potential negative overall\n",
    "    * installed capacity\n",
    "    * fixed and variable costs\n",
    "    * specific investments"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Introduce data fixes\n",
    "Steps applied:\n",
    "* Choose potential and cost parameters to further analyse\n",
    "* Perform **interpolation on potential and capital costs parameters** in order to\n",
    "    1. Fill data gaps and\n",
    "    2. remove inplausibilities such as potentials changing very strong and not consistent within the five year intervalls used.\n",
    "* Combine the data to an overall data set once the actions described in the following are done.\n",
    "\n",
    "In order to come up with consistent trends, it is proceeded as follows:\n",
    "* For *potential-related* parameters, values for the status quo, 2030 and 2050 are assessed. This serves to\n",
    "    * depict intermediate trends (such as a temporary increase of a technology) and\n",
    "    * create a consistent long-term development.\n",
    "* For *potential-related* parameters, a strong interlinkage exists. In order not to create inconsistencies,\n",
    "    * The median ratio between negative and positive shifting potentials is determined and kept constant over time.\n",
    "    * The same holds for the median ratio between shedding and installed capacity.\n",
    "    * Before, it is ensured that installed capacity is always the maximum value of the capacity-related parameters.\n",
    "    * These corrections serve to prevent data inconsistensies over time, such as a strongly increasing positive shifting potential while installed capacity only weakly increases or even decreases.\n",
    "    * It is worth noting that this ultimately affects potential results. Thus, the trade-off between keeping an inconsistent data basis over time and having a consistent one, but with altered technical potential values has been decided in favor of consistency. The author wants to stress that the high uncertainties for potentials in combination with contradictory literature as well as the intention to keep as much of the original data basis as possible and to shed light on the uncertainty range (using 5%, 50% and 95% percentile values) serve as a justification for the chosen approach.\n",
    "* For *investment expenses* parameters, the minimum cost value given is assigned to 2050, thus assuming a cost reduction over time. Zero cost values are replaced by 10% of the highest value to depict cost degression over time. Also, a minimum value of 0.01 €/kW is assigned.\n",
    "* For *fixed costs*, these are assumed as a percentage share of investment expenses. 2% of investment expenses or at least 0.01 €/kW * a are used. A variation is included since the investment expenses vary by demand response scenario.\n",
    "* For *variable costs*, use mean values across all years and assume costs to remain constant (in real terms)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# choose parameters to be used (parameters for which some adaptions are needed)\n",
    "params_to_use = [\n",
    "    \"potential_neg_overall\",\n",
    "    \"potential_pos_overall\",\n",
    "    \"potential_pos_overall_shed\",\n",
    "    \"installed_cap\",\n",
    "    \"fixed_costs\",\n",
    "    \"variable_costs\",\n",
    "    \"variable_costs_shed\",\n",
    "    \"specific_investments\"\n",
    "]\n",
    "\n",
    "params_to_use = [(a, b) for a in quantile_cols for b in params_to_use]\n",
    "\n",
    "# Slice the parameter values needed\n",
    "slice_params = all_params_reshaped.loc[:, params_to_use]\n",
    "slice_params_unjoined = all_params_reshaped_unjoined.loc[:, params_to_use]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define potential and cost cols\n",
    "pot_cols = [\n",
    "    \"potential_neg_overall\",\n",
    "    \"potential_pos_overall\", \n",
    "    \"potential_pos_overall_shed\",\n",
    "    \"installed_cap\"\n",
    "]\n",
    "pot_cols = [(a, b) for a in quantile_cols for b in pot_cols]\n",
    "\n",
    "invest_cols = [\n",
    "    \"specific_investments\",\n",
    "\n",
    "]\n",
    "invest_cols = [(a, b) for a in quantile_cols for b in invest_cols]\n",
    "\n",
    "fixed_costs_cols = [\n",
    "    \"fixed_costs\"\n",
    "]\n",
    "fixed_costs_cols = [(a, b) for a in quantile_cols for b in fixed_costs_cols]\n",
    "\n",
    "variable_costs_cols = [\n",
    "    \"variable_costs\",\n",
    "    \"variable_costs_shed\"   \n",
    "]\n",
    "variable_costs_cols = [(a, b) for a in quantile_cols for b in variable_costs_cols]\n",
    "\n",
    "# Determine demand response categories to use for assigning values:\n",
    "# Use first two index levels corresponding to process category and sector\n",
    "dr_categories = all_params_reshaped.reset_index(level=2).index.unique()\n",
    "dr_categories_unjoined = all_params_reshaped_unjoined.reset_index(level=2).index.unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Combined data set\n",
    "# Instanciate new DataFrame to store manipulated outputs\n",
    "parameters_for_clustering = pd.DataFrame(\n",
    "    index=pd.MultiIndex.from_product([[], [], [], []],\n",
    "    names=[\"Prozesskategorie\", \"sector\", \"year\", \"kind\"]),\n",
    "    columns=pd.MultiIndex.from_product([[], []])\n",
    ")\n",
    "\n",
    "# Create a list of DataFrames to concat\n",
    "to_concat = [parameters_for_clustering]\n",
    "\n",
    "original_pot_cols = [\n",
    "    \"potential_neg_overall\",\n",
    "    \"potential_pos_overall\", \n",
    "    \"potential_pos_overall_shed\",\n",
    "    \"installed_cap\"\n",
    "]\n",
    "\n",
    "for category in dr_categories:\n",
    "    process = category[0]\n",
    "    sector = category[1]\n",
    "    kind = category[2]\n",
    "    \n",
    "    multi_ix = pd.MultiIndex.from_product(\n",
    "        [[process], [sector], years, [kind]], \n",
    "        names=[\"Prozesskategorie\", \"sector\", \"year\", \"kind\"]\n",
    "    )\n",
    "\n",
    "    new_df = pd.DataFrame(\n",
    "        index=multi_ix, \n",
    "        columns=pd.MultiIndex.from_tuples(pot_cols + invest_cols + fixed_costs_cols + variable_costs_cols)\n",
    "    )\n",
    "    \n",
    "    for quantile in quantile_cols:\n",
    "        cols_to_use = {\n",
    "            \"pot\": [(quantile, col) for col in original_pot_cols],\n",
    "            \"pot_pos\": [(quantile, \"potential_pos_overall\")],\n",
    "            \"pot_neg\": [(quantile, \"potential_neg_overall\")],\n",
    "            \"pot_shed\": [(quantile, \"potential_pos_overall_shed\")],\n",
    "            \"installed_cap\": [(quantile, \"installed_cap\")]\n",
    "        }\n",
    "        potential_values = {\"SQ\": {}, \"2030\": {}, \"2050\": {}}\n",
    "\n",
    "        # Use potential values for status quo, 2030 (if available) and 2050\n",
    "        try:\n",
    "            potential_values[\"SQ\"] = {\n",
    "                key: slice_params.loc[(process, sector, \"SQ\", kind), col].values \n",
    "                for key, col in cols_to_use.items()\n",
    "            }\n",
    "        except KeyError:\n",
    "            continue\n",
    "        try:\n",
    "            potential_values[\"2030\"] = {\n",
    "                key: slice_params.loc[(process, sector, \"2030\", kind), col].values \n",
    "                for key, col in cols_to_use.items()\n",
    "            }\n",
    "        except:\n",
    "            array = np.empty((1,))\n",
    "            array.fill(np.nan)\n",
    "            potential_values[\"2030\"] = {\n",
    "                key: array for key in cols_to_use\n",
    "            }\n",
    "        try:\n",
    "            potential_values[\"2050\"] = {\n",
    "                key: slice_params.loc[(process, sector, \"2050\", kind), col].values \n",
    "                for key, col in cols_to_use.items()\n",
    "            }\n",
    "        except:\n",
    "            potential_values[\"2050\"] = potential_values[\"SQ\"]\n",
    "            \n",
    "        # Determine mean ratios between negative and positive potential\n",
    "        all_ratios_neg_to_pos = np.concatenate(\n",
    "            [potential_values[y][\"pot_neg\"] / potential_values[y][\"pot_pos\"] for y in potential_values]\n",
    "        )\n",
    "        all_ratios_shed_to_installed = np.concatenate(\n",
    "            [potential_values[y][\"pot_shed\"] / potential_values[y][\"installed_cap\"] for y in potential_values]\n",
    "        )\n",
    "        # Replace infinity values (division by zero)\n",
    "        all_ratios_neg_to_pos[all_ratios_neg_to_pos == np.inf] = np.nan\n",
    "        all_ratios_shed_to_installed[all_ratios_shed_to_installed == np.inf] = np.nan\n",
    "        ratio_neg_to_pos_potential_data = np.nanmean(all_ratios_neg_to_pos)\n",
    "        ratio_shed_to_installed_data = np.nanmean(all_ratios_shed_to_installed)  \n",
    "\n",
    "        # Set potential values\n",
    "        for col, col_name in cols_to_use.items():\n",
    "            if col != \"pot\":\n",
    "                for y in potential_values:\n",
    "                    new_df.loc[(process, sector, y, kind), col_name] = potential_values[y][col]\n",
    "        \n",
    "        # Do some corrections\n",
    "        for y in potential_values:\n",
    "            # Use ratio to adjust negative potential values using mean ratio between negative and positive potential\n",
    "            new_df.loc[(process, sector, y, kind), cols_to_use[\"pot_neg\"]] = (\n",
    "                new_df.loc[(process, sector, y, kind), cols_to_use[\"pot_pos\"]] * ratio_neg_to_pos_potential_data\n",
    "            ).values\n",
    "            # Ensure installed capacity to be the highest value\n",
    "            new_df.loc[(process, sector, y, kind), cols_to_use[\"installed_cap\"]] = np.max(\n",
    "                new_df.loc[(process, sector, y, kind), cols_to_use[\"pot\"]]\n",
    "            )\n",
    "            # Use ratio to adjust shedding potential values using mean ratio between \n",
    "            # shedding potential and installed capacity  \n",
    "            new_df.loc[(process, sector, y, kind), cols_to_use[\"pot_shed\"]] = (\n",
    "                new_df.loc[(process, sector, y, kind), cols_to_use[\"installed_cap\"]] * ratio_shed_to_installed_data\n",
    "            ).values\n",
    "    \n",
    "    # For variable costs values, use mean values across all years, ignoring duplicate values\n",
    "    variable_costs_vals = slice_params.loc[(process, sector, years, kind), variable_costs_cols]\n",
    "    for col in variable_costs_vals.columns:\n",
    "        variable_costs_vals[col] = np.mean(variable_costs_vals[col].unique())  \n",
    "    new_df.loc[(process, sector, years, kind), variable_costs_cols] = variable_costs_vals.replace({0: 0.01})\n",
    "    \n",
    "    # Use investment expenses values\n",
    "    invest_vals_SQ = slice_params.loc[(process, sector, \"SQ\", kind), invest_cols]\n",
    "    invest_vals_SQ = invest_vals_SQ.replace({0: 0.01})\n",
    "    min_costs = slice_params.loc[(process, sector, years, kind), invest_cols].min()\n",
    "    min_costs = np.maximum(min_costs, 0.1 * invest_vals_SQ)\n",
    "    invest_vals_2050 = np.where(min_costs < 0.01, 0.01, min_costs)    \n",
    "    \n",
    "    # Assign minimum investment expenses value to the year 2050\n",
    "    # Assume investment expenses not to fall below 10% of original value\n",
    "    new_df.loc[(process, sector, \"SQ\", kind), invest_cols] = invest_vals_SQ.fillna(0.01)\n",
    "    new_df.loc[(process, sector, \"2050\", kind), invest_cols] = invest_vals_2050\n",
    "    \n",
    "    # Set fixed costs to 2% of investment expenses or at least 0.01 €/kW * a\n",
    "    for quantile in quantile_cols:\n",
    "        for y in [\"SQ\", \"2050\"]:\n",
    "            fixed_costs_val = new_df.loc[\n",
    "                (process, sector, y, kind), [(quantile, \"specific_investments\")]\n",
    "            ].values * 0.02\n",
    "            new_df.loc[(process, sector, y, kind), [(quantile, \"fixed_costs\")]] = np.where(\n",
    "                fixed_costs_val < 0.01, 0.01, fixed_costs_val\n",
    "            )\n",
    "    \n",
    "    # Correct dtype, interpolate between status quo, 2030 and 2050 and create common data basis again\n",
    "    new_df = new_df.astype(\"float64\")\n",
    "    new_df = new_df.interpolate(axis=0)\n",
    "    \n",
    "    to_concat.append(new_df)\n",
    "\n",
    "parameters_for_clustering = pd.concat([el for el in to_concat], levels=([0, 1, 2]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Uncombined data set\n",
    "# Instanciate new DataFrame to store manipulated outputs\n",
    "parameters_for_clustering_unjoined = pd.DataFrame(\n",
    "    index=pd.MultiIndex.from_product([[], [], [], []],\n",
    "    names=[\"Prozesskategorie\", \"sector\", \"year\", \"kind\"]),\n",
    "    columns=pd.MultiIndex.from_product([[], []])\n",
    ")\n",
    "\n",
    "# Create a list of DataFrames to concat\n",
    "to_concat = [parameters_for_clustering_unjoined]\n",
    "\n",
    "original_pot_cols = [\n",
    "    \"potential_neg_overall\",\n",
    "    \"potential_pos_overall\", \n",
    "    \"potential_pos_overall_shed\",\n",
    "    \"installed_cap\"\n",
    "]\n",
    "\n",
    "for category in dr_categories_unjoined:\n",
    "    process = category[0]\n",
    "    sector = category[1]\n",
    "    kind = category[2]\n",
    "    \n",
    "    multi_ix = pd.MultiIndex.from_product(\n",
    "        [[process], [sector], years, [kind]], \n",
    "        names=[\"Prozesskategorie\", \"sector\", \"year\", \"kind\"]\n",
    "    )\n",
    "\n",
    "    new_df = pd.DataFrame(\n",
    "        index=multi_ix, \n",
    "        columns=pd.MultiIndex.from_tuples(pot_cols + invest_cols + fixed_costs_cols + variable_costs_cols)\n",
    "    )\n",
    "    \n",
    "    for quantile in quantile_cols:\n",
    "        cols_to_use = {\n",
    "            \"pot\": [(quantile, col) for col in original_pot_cols],\n",
    "            \"pot_pos\": [(quantile, \"potential_pos_overall\")],\n",
    "            \"pot_neg\": [(quantile, \"potential_neg_overall\")],\n",
    "            \"pot_shed\": [(quantile, \"potential_pos_overall_shed\")],\n",
    "            \"installed_cap\": [(quantile, \"installed_cap\")]\n",
    "        }\n",
    "        potential_values = {\"SQ\": {}, \"2030\": {}, \"2050\": {}}\n",
    "        \n",
    "        # Use potential values for status quo, 2030 (if available) and 2050\n",
    "        try:\n",
    "            potential_values[\"SQ\"] = {\n",
    "                key: slice_params_unjoined.loc[(process, sector, \"SQ\", kind), col].values \n",
    "                for key, col in cols_to_use.items()\n",
    "            }\n",
    "        except KeyError:\n",
    "            continue\n",
    "        try:\n",
    "            potential_values[\"2030\"] = {\n",
    "                key: slice_params_unjoined.loc[(process, sector, \"2030\", kind), col].values \n",
    "                for key, col in cols_to_use.items()\n",
    "            }\n",
    "        except:\n",
    "            array = np.empty((1,))\n",
    "            array.fill(np.nan)\n",
    "            potential_values[\"2030\"] = {\n",
    "                key: array for key in cols_to_use\n",
    "            }\n",
    "        try:\n",
    "            potential_values[\"2050\"] = {\n",
    "                key: slice_params_unjoined.loc[(process, sector, \"2050\", kind), col].values \n",
    "                for key, col in cols_to_use.items()\n",
    "            }\n",
    "        except:\n",
    "            potential_values[\"2050\"] = potential_values[\"SQ\"]\n",
    "\n",
    "        # Determine mean ratios between negative and positive potential\n",
    "        all_ratios_neg_to_pos = np.concatenate(\n",
    "            [potential_values[y][\"pot_neg\"] / potential_values[y][\"pot_pos\"] for y in potential_values]\n",
    "        )\n",
    "        all_ratios_shed_to_installed = np.concatenate(\n",
    "            [potential_values[y][\"pot_shed\"] / potential_values[y][\"installed_cap\"] for y in potential_values]\n",
    "        )\n",
    "        # Replace infinity values (division by zero)\n",
    "        all_ratios_neg_to_pos[all_ratios_neg_to_pos == np.inf] = np.nan\n",
    "        all_ratios_shed_to_installed[all_ratios_shed_to_installed == np.inf] = np.nan\n",
    "        ratio_neg_to_pos_potential_data = np.nanmean(all_ratios_neg_to_pos)\n",
    "        ratio_shed_to_installed_data = np.nanmean(all_ratios_shed_to_installed)  \n",
    "\n",
    "        # Set potential values\n",
    "        for col, col_name in cols_to_use.items():\n",
    "            if col != \"pot\":\n",
    "                for y in potential_values:\n",
    "                    new_df.loc[(process, sector, y, kind), col_name] = potential_values[y][col]\n",
    "        \n",
    "        # Do some corrections\n",
    "        for y in potential_values:\n",
    "            # Use ratio to adjust negative potential values using mean ratio between negative and positive potential\n",
    "            new_df.loc[(process, sector, y, kind), cols_to_use[\"pot_neg\"]] = (\n",
    "                new_df.loc[(process, sector, y, kind), cols_to_use[\"pot_pos\"]] * ratio_neg_to_pos_potential_data\n",
    "            ).values\n",
    "            # Ensure installed capacity to be the highest value\n",
    "            new_df.loc[(process, sector, y, kind), cols_to_use[\"installed_cap\"]] = np.max(\n",
    "                new_df.loc[(process, sector, y, kind), cols_to_use[\"pot\"]]\n",
    "            )\n",
    "            # Use ratio to adjust shedding potential values using mean ratio between \n",
    "            # shedding potential and installed capacity  \n",
    "            new_df.loc[(process, sector, y, kind), cols_to_use[\"pot_shed\"]] = (\n",
    "                new_df.loc[(process, sector, y, kind), cols_to_use[\"installed_cap\"]] * ratio_shed_to_installed_data\n",
    "            ).values\n",
    "\n",
    "    # For variable costs values, use mean values across all years, ignoring duplicate values\n",
    "    variable_costs_vals = slice_params_unjoined.loc[(process, sector, years, kind), variable_costs_cols]\n",
    "    for col in variable_costs_vals.columns:\n",
    "        variable_costs_vals[col] = np.mean(variable_costs_vals[col].unique())  \n",
    "    new_df.loc[(process, sector, years, kind), variable_costs_cols] = variable_costs_vals.replace({0: 0.01})\n",
    "    \n",
    "    # Use investment expenses values\n",
    "    invest_vals_SQ = slice_params_unjoined.loc[(process, sector, \"SQ\", kind), invest_cols]\n",
    "    invest_vals_SQ = invest_vals_SQ.replace({0: 0.01})\n",
    "    min_costs = slice_params_unjoined.loc[(process, sector, years, kind), invest_cols].min()\n",
    "    min_costs = np.maximum(min_costs, 0.1 * invest_vals_SQ)\n",
    "    invest_vals_2050 = np.where(min_costs < 0.01, 0.01, min_costs)\n",
    "\n",
    "    # Assign minimum investment expenses value to the year 2050\n",
    "    # Assume investment expenses not to fall below 10% of original value\n",
    "    new_df.loc[(process, sector, \"SQ\", kind), invest_cols] = invest_vals_SQ.fillna(0.01)\n",
    "    new_df.loc[(process, sector, \"2050\", kind), invest_cols] = invest_vals_2050\n",
    "    \n",
    "    # Set fixed costs to 2% of investment expenses or at least 0.01 €/kW * a\n",
    "    for quantile in quantile_cols:\n",
    "        for y in [\"SQ\", \"2050\"]:\n",
    "            fixed_costs_val = new_df.loc[\n",
    "                (process, sector, y, kind), [(quantile, \"specific_investments\")]\n",
    "            ].values * 0.02\n",
    "            new_df.loc[(process, sector, y, kind), [(quantile, \"fixed_costs\")]] = np.where(\n",
    "                fixed_costs_val < 0.01, 0.01, fixed_costs_val\n",
    "            )\n",
    "    \n",
    "    # Correct dtype, interpolate between status quo, 2030 and 2050 and create common data basis again\n",
    "    new_df = new_df.astype(\"float64\")\n",
    "    new_df = new_df.interpolate(axis=0)\n",
    "    \n",
    "    to_concat.append(new_df)\n",
    "\n",
    "parameters_for_clustering_unjoined = pd.concat([el for el in to_concat], levels=([0, 1, 2]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Assignment of remaining (time-related parameters)\n",
    "Assume time-related parameters to remain constant over time, i.e. the same as in the status quo data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Assign remaining (mostly time-related) parameters\n",
    "params_remaining = [\n",
    "    \"interference_duration_neg\",\n",
    "    \"interference_duration_pos\",\n",
    "    \"interference_duration_pos_shed\",\n",
    "    \"maximum_activations_year\",\n",
    "    \"maximum_activations_year_shed\",\n",
    "    \"regeneration_duration\", \n",
    "    \"shifting_duration\"\n",
    "]\n",
    "\n",
    "params_remaining = [\n",
    "    (a, b) for a in quantile_cols for b in params_remaining\n",
    "]\n",
    "\n",
    "# Slice the parameter values needed (both data sets - combined and uncombined)\n",
    "slice_params = all_params_reshaped.loc[:, params_remaining]\n",
    "slice_params_unjoined = all_params_reshaped_unjoined.loc[:, params_remaining]\n",
    "\n",
    "parameters_for_clustering = parameters_for_clustering.reindex(\n",
    "    columns=pot_cols + invest_cols + fixed_costs_cols + variable_costs_cols + params_remaining\n",
    ")\n",
    "parameters_for_clustering_unjoined = parameters_for_clustering_unjoined.reindex(\n",
    "    columns=pot_cols + invest_cols + fixed_costs_cols + variable_costs_cols + params_remaining\n",
    ")\n",
    "\n",
    "# Use status quo values (best data basis) for all years\n",
    "parameters_for_clustering.loc[\n",
    "    parameters_for_clustering.index.get_level_values(2) == \"SQ\", params_remaining\n",
    "] = (\n",
    "    slice_params.loc[slice_params.index.get_level_values(2) == \"SQ\", params_remaining]\n",
    ")\n",
    "parameters_for_clustering.fillna(method=\"ffill\", inplace=True)\n",
    "\n",
    "parameters_for_clustering_unjoined.loc[\n",
    "    parameters_for_clustering_unjoined.index.get_level_values(2) == \"SQ\", params_remaining\n",
    "] = (\n",
    "    slice_params_unjoined.loc[slice_params_unjoined.index.get_level_values(2)\n",
    "                              == \"SQ\", params_remaining]\n",
    ")\n",
    "parameters_for_clustering_unjoined.fillna(method=\"ffill\", inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Slightly rearrange and filter status quo data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Flatten columns (both data sets - combined and uncombined)\n",
    "parameters_for_clustering.loc[(\"new_cols\",\"\",\"\",\"\"), :] = (\n",
    "    parameters_for_clustering.columns.get_level_values(0) \n",
    "    + \"_\" + parameters_for_clustering.columns.get_level_values(1)\n",
    ")\n",
    "\n",
    "parameters_for_clustering.columns = parameters_for_clustering.loc[(\"new_cols\",\"\",\"\",\"\")]\n",
    "parameters_for_clustering.columns.name = \"parameters\"\n",
    "parameters_for_clustering.drop(index=(\"new_cols\",\"\",\"\",\"\"), inplace=True)\n",
    "\n",
    "parameters_for_clustering_unjoined.loc[(\"new_cols\",\"\",\"\",\"\"), :] = (\n",
    "    parameters_for_clustering_unjoined.columns.get_level_values(0) \n",
    "    + \"_\" + parameters_for_clustering_unjoined.columns.get_level_values(1)\n",
    ")\n",
    "\n",
    "parameters_for_clustering_unjoined.columns = parameters_for_clustering_unjoined.loc[(\"new_cols\",\"\",\"\",\"\")]\n",
    "parameters_for_clustering_unjoined.columns.name = \"parameters\"\n",
    "parameters_for_clustering_unjoined.drop(index=(\"new_cols\",\"\",\"\",\"\"), inplace=True)\n",
    "\n",
    "# Filter status quo data\n",
    "parameters_for_clustering_status_quo = parameters_for_clustering.loc[\n",
    "    parameters_for_clustering.index.get_level_values(2) == \"SQ\"\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cluster data\n",
    "* Do a clustering using K-Means and the cluster parameters defined in the above parameter settings (alternative: agglomerative clustering using ward linkage).\n",
    "* Default: Cluster within sectors by shifting times, positive interference duration, variable costs as well as specific investments"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Do the actual clustering\n",
    "* Determine the number of clusters per sector and shifting resp. shedding eligibility\n",
    "* Perform a K-Means clustering on the data (alternative: agglomerative clustering using ward linkage)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize the column to store the cluster label\n",
    "parameters_for_clustering_status_quo[\"cluster\"] = 0\n",
    "\n",
    "cluster_string = f\"Number of clusters for {cluster_algo} clustering:\"\n",
    "print(cluster_string)\n",
    "print(\"-\" * len(cluster_string)) \n",
    "\n",
    "# Increase cluster numbers by (arbitratily chosen) increments\n",
    "# in order not to overwrite cluster information\n",
    "increment = 100\n",
    "\n",
    "for sector, kind in parameters_for_clustering_status_quo.reset_index(level=[0, 2]).index.unique():\n",
    "    # Determine the number of clusters per sector and print it\n",
    "    n = math.ceil(\n",
    "        parameters_for_clustering_status_quo[\n",
    "            (parameters_for_clustering_status_quo.index.get_level_values(1) == sector)\n",
    "            & (parameters_for_clustering_status_quo.index.get_level_values(3) == kind)\n",
    "        ].shape[0] * share_clusters\n",
    "    )\n",
    "    print(f\"{sector: <9} / {kind : <20}: {n : >3}\")\n",
    "    \n",
    "    if cluster_algo == \"KMeans\":\n",
    "        # Do the actual clustering and assign the cluster labels\n",
    "        parameters_for_clustering_status_quo.loc[\n",
    "            (parameters_for_clustering_status_quo.index.get_level_values(1) == sector)\n",
    "            & (parameters_for_clustering_status_quo.index.get_level_values(3) == kind), \n",
    "            \"cluster\"] = (\n",
    "                KMeans(n_clusters=n).fit(\n",
    "                    parameters_for_clustering_status_quo.loc[\n",
    "                        (parameters_for_clustering_status_quo.index.get_level_values(1) == sector)\n",
    "                        & (parameters_for_clustering_status_quo.index.get_level_values(3) == kind), \n",
    "                        cluster_parameters\n",
    "                    ].values\n",
    "                ).labels_ + 1\n",
    "            )\n",
    "        \n",
    "    elif cluster_algo == \"ward\":\n",
    "        # Do the actual clustering and assign the cluster labels\n",
    "        parameters_for_clustering_status_quo.loc[\n",
    "            parameters_for_clustering_status_quo.index.get_level_values(1) == sector, \"cluster\"] = (\n",
    "                AgglomerativeClustering(n_clusters=n, linkage=\"ward\").fit(\n",
    "                    parameters_for_clustering_status_quo.loc[\n",
    "                        parameters_for_clustering_status_quo.index.get_level_values(1) \n",
    "                        == sector, cluster_parameters\n",
    "                    ].values\n",
    "                ).labels_ + 1       \n",
    "            )\n",
    "\n",
    "    # Increment the cluster labels to avoid overwriting in the next iteration \n",
    "    # (cluster numbers start with zero)\n",
    "    parameters_for_clustering_status_quo.loc[\n",
    "        parameters_for_clustering_status_quo[\"cluster\"] != 0, \"cluster\"\n",
    "    ] += increment\n",
    "\n",
    "# Print the number of unique cluster labels (cross check)\n",
    "print(f\"Original number of clusters     : {parameters_for_clustering_status_quo['cluster'].nunique() : >3}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Cluster assignment**:\n",
    "Assign cluster information for future years"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Drop year from index\n",
    "parameters_for_clustering.reset_index(level=2, inplace=True)\n",
    "parameters_for_clustering_unjoined.reset_index(level=2, inplace=True)\n",
    "parameters_for_clustering_status_quo.reset_index(level=2, inplace=True)\n",
    "\n",
    "for year in years:\n",
    "    parameters_for_clustering.loc[\n",
    "        parameters_for_clustering[\"year\"] == year, \"cluster\"\n",
    "    ] = parameters_for_clustering_status_quo[\"cluster\"].values\n",
    "    \n",
    "    # Assgin values for uncombined data set\n",
    "    parameters_for_clustering_unjoined.loc[\n",
    "        parameters_for_clustering_unjoined[\"year\"] == year, \"cluster\"\n",
    "    ] = parameters_for_clustering_status_quo[\"cluster\"]\n",
    "    \n",
    "    # Handle missing cluster information\n",
    "    no_cluster_ixs = list(set(parameters_for_clustering_unjoined[\n",
    "        (parameters_for_clustering_unjoined[\"year\"] == year)\n",
    "        & (parameters_for_clustering_unjoined[\"cluster\"].isna())\n",
    "    ].index.get_level_values(0)))\n",
    "    no_cluster_ixs_tcs_hoho = [(i, j, k) for i in no_cluster_ixs for j in [\"tcs+hoho\"] for k in [\"shift_only\"]]\n",
    "\n",
    "    # Use only the first index level to be able to assign values (second index levels won't match)\n",
    "    clusters_to_use = parameters_for_clustering[\n",
    "        parameters_for_clustering[\"year\"] == year \n",
    "    ].loc[no_cluster_ixs_tcs_hoho, \"cluster\"].to_frame()\n",
    "    clusters_to_use = clusters_to_use.set_index(clusters_to_use.index.get_level_values(0))\n",
    "\n",
    "    # Use an excerpt of the overall DataFrame to assign the cluster info\n",
    "    excerpt = parameters_for_clustering_unjoined.loc[\n",
    "        (parameters_for_clustering_unjoined[\"year\"] == year)\n",
    "        & (parameters_for_clustering_unjoined[\"cluster\"].isna()),\n",
    "        \"cluster\"\n",
    "    ].to_frame()\n",
    "    excerpt = excerpt.reset_index().set_index(\"Prozesskategorie\")\n",
    "\n",
    "    excerpt[\"cluster\"] = clusters_to_use[\"cluster\"]\n",
    "    excerpt = excerpt.set_index(\"sector\", append=True)\n",
    "\n",
    "    parameters_for_clustering_unjoined.loc[\n",
    "        (parameters_for_clustering_unjoined[\"year\"] == year)\n",
    "        & (parameters_for_clustering_unjoined[\"cluster\"].isna()),\n",
    "        \"help_sector\"\n",
    "    ] = \"tcs+hoho\"\n",
    "    parameters_for_clustering_unjoined.loc[\n",
    "        (parameters_for_clustering_unjoined[\"year\"] == year)\n",
    "        & (parameters_for_clustering_unjoined[\"cluster\"].isna()),\n",
    "        \"cluster\"\n",
    "    ] = excerpt.loc[:, \"cluster\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Show and explore the clusters\n",
    "Print the clusters in order to visually inspect them"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": false
    }
   },
   "outputs": [],
   "source": [
    "if print_clusters:\n",
    "    for el in np.sort(parameters_for_clustering_status_quo[\"cluster\"].unique()):\n",
    "        print(20 * \"-\")\n",
    "        print(\"cluster number: \"+str(el))\n",
    "        print(20 * \"-\")\n",
    "        display(parameters_for_clustering_status_quo[parameters_for_clustering_status_quo[\"cluster\"] == el])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Drop unsensible tcs cluster:\n",
    "* Lighting in the tcs sector (flower industry / greenhouses) is supposed to be eligible for load shedding, but does not have any shedding potential information.\n",
    "* Thus, the process and its corresponding cluster will be removed from the data set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Drop kind (shifting / shedding) from index\n",
    "parameters_for_clustering = parameters_for_clustering.reset_index(level=2, drop=False)\n",
    "parameters_for_clustering_unjoined = parameters_for_clustering_unjoined.reset_index(level=2, drop=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Drop unsensible tcs cluster\n",
    "to_drop = [(\"Beleuchtung\", \"tcs\")]\n",
    "parameters_for_clustering_status_quo.drop(index=to_drop, inplace=True)\n",
    "parameters_for_clustering.drop(index=to_drop, inplace=True)\n",
    "parameters_for_clustering_unjoined.drop(index=to_drop, inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "is_executing": false,
     "name": "#%% md\n"
    },
    "toc-hr-collapsed": false
   },
   "source": [
    "## Read in and match availability data\n",
    "* Read in the availability time series data collected from Benz (2019), Odeh (2019) and Stange (2019).\n",
    "* Match the data to the respective demand response categories."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "is_executing": false,
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# Read in availability time series data\n",
    "availability_ind_pos = pd.read_excel(path_folder_availability+filename_availability_in,\n",
    "                                     sheet_name = \"ind_pos\", header=0, index_col=0)\n",
    "\n",
    "availability_ind_neg = pd.read_excel(path_folder_availability+filename_availability_in,\n",
    "                                     sheet_name = \"ind_neg\", header=0, index_col=0)\n",
    "\n",
    "availability_tcs_pos = pd.read_excel(path_folder_availability+filename_availability_in,\n",
    "                                     sheet_name = \"tcs_pos\", header=0, index_col=0)\n",
    "\n",
    "availability_tcs_neg = pd.read_excel(path_folder_availability+filename_availability_in,\n",
    "                                     sheet_name = \"tcs_neg\", header=0, index_col=0)\n",
    "\n",
    "availability_hoho_pos = pd.read_excel(path_folder_availability+filename_availability_in,\n",
    "                                     sheet_name = \"hoho_pos\", header=0, index_col=0)\n",
    "\n",
    "availability_hoho_neg = pd.read_excel(path_folder_availability+filename_availability_in,\n",
    "                                     sheet_name = \"hoho_neg\", header=0, index_col=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "is_executing": false,
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# Map the demand response categories to the column names used for the availability time series\n",
    "# (naming conventions from the bachelor theses are widely used)\n",
    "availability_cats_ind_pos = {\n",
    "    'normiertes\\nLRP Alu': ('Primäraluminiumelektrolyse', 'ind'),\n",
    "    'normiertes\\nLRP CHL': ('Chlor-Alkali-Elektrolyse', 'ind'),\n",
    "    'normiertes\\nLRP HS': ('Holz- und Zellstoffherstellung', 'ind'),\n",
    "    'normiertes\\nLRP PM': ('Papiermaschinen', 'ind'),\n",
    "    'normiertes\\nLRP ES': ('Elektrostahlherstellung (Lichtbogenofen)', 'ind'),\n",
    "    'normiertes\\nLRP RZM': ('Zementherstellung', 'ind'),\n",
    "    'normiertes\\nLRP KUZI': ('Kupfer- und Zinkherstellung (Elektrolyse)', 'ind'),\n",
    "    'normiertes\\nLRP AP': ('Altpapierrecycling (Pulper)', 'ind')\n",
    "}\n",
    "\n",
    "availability_cats_ind_neg = {\n",
    "    'normiertes\\nLZP Alu': ('Primäraluminiumelektrolyse', 'ind'),\n",
    "    'normiertes\\nLZP CHL': ('Chlor-Alkali-Elektrolyse', 'ind'),\n",
    "    'normiertes\\nLZP HS': ('Holz- und Zellstoffherstellung', 'ind'),\n",
    "    'normiertes\\nLZP PM': ('Papiermaschinen', 'ind'),\n",
    "    'normiertes\\nLZP ES': ('Elektrostahlherstellung (Lichtbogenofen)', 'ind'),\n",
    "    'normiertes\\nLZP RZM': ('Zementherstellung', 'ind'),\n",
    "    'normiertes\\nLZP KUZI': ('Kupfer- und Zinkherstellung (Elektrolyse)', 'ind'),\n",
    "    'normiertes\\nLZP AP': ('Altpapierrecycling (Pulper)', 'ind')\n",
    "}\n",
    "\n",
    "# Note: Instead of KGR in the tcs sector, profiles for food retailing are introduced and used.\n",
    "availability_cats_tcs_pos = {\n",
    "    'Lastabschaltung LÜ normiert': ('Belüftung', 'tcs'), \n",
    "    'Lastabschaltung WVP normiert': ('Pumpenanwendungen in der Wasserversorgung', 'tcs'),\n",
    "    'Lastabschaltung KGR normiert': ('Prozesskälte Handel', 'tcs'),\n",
    "    'Lastabschaltung KÜ normiert': ('Kühlhäuser', 'tcs'),\n",
    "    'Lastabschaltung KA normiert': ('Klimakälte', 'tcs'),  \n",
    "    'Lastabschaltung EH normiert': ('Nachtspeicherheizungen', 'tcs'),\n",
    "    'Lastabschaltung WP normiert': ('Wärmepumpen', 'tcs'), \n",
    "    'Lastabschaltung WW normiert': ('Warmwasserbereitstellung', 'tcs')\n",
    "}\n",
    "\n",
    "availability_cats_tcs_neg = {\n",
    "    'Lastzuschaltung LÜ normiert': ('Belüftung', 'tcs'),\n",
    "    'Lastzuschaltung WVP normiert': ('Pumpenanwendungen in der Wasserversorgung', 'tcs'),\n",
    "    'Lastzuschaltung KGR normiert': ('Prozesskälte Handel', 'tcs'),\n",
    "    'Lastzuschaltung KÜ normiert': ('Kühlhäuser', 'tcs'),\n",
    "    'Lastzuschaltung KA normiert': ('Klimakälte', 'tcs'), \n",
    "    'Lastzuschaltung EH normiert': ('Nachtspeicherheizungen', 'tcs'),\n",
    "    'Lastzuschaltung WP normiert': ('Wärmepumpen', 'tcs'), \n",
    "    'Lastzuschaltung WW normiert': ('Warmwasserbereitstellung', 'tcs')\n",
    "}\n",
    "\n",
    "availability_cats_hoho_pos = {\n",
    "    'Lastabschaltung KGR': ('Kühlschränke', 'hoho'), \n",
    "    'Lastabschaltung WM': ('Waschmaschinen', 'hoho'), \n",
    "    'Lastabschaltung WT': ('Wäschetrockner', 'hoho'),\n",
    "    'Lastabschaltung GS': ('Geschirrspüler', 'hoho'), \n",
    "    'Lastabschaltung NH normiert': ('Nachtspeicherheizungen', 'hoho'),\n",
    "    'Lastabschaltung WP normiert': ('Wärmepumpen', 'hoho'), \n",
    "    'Lastabschaltung UP': ('Heizungsumwälzpumpen', 'hoho'),\n",
    "    'Lastabschaltung RK': ('Klimakälte', 'hoho'), \n",
    "    'Lastabschaltung WW Tag': ('Warmwasserbereitstellung', 'hoho')\n",
    "}\n",
    "\n",
    "availability_cats_hoho_neg = {\n",
    "    'Lastzuschaltung KGR normiert': ('Kühlschränke', 'hoho'),\n",
    "    'Lastzuschaltung WM normiert': ('Waschmaschinen', 'hoho'),\n",
    "    'Lastzuschaltung WT normiert': ('Wäschetrockner', 'hoho'),\n",
    "    'Lastzuschaltung GS normiert': ('Geschirrspüler', 'hoho'), \n",
    "    'Lastzuschaltung NH normiert': ('Nachtspeicherheizungen', 'hoho'),\n",
    "    'Lastzuschaltung WP normiert': ('Wärmepumpen', 'hoho'),\n",
    "    'Lastzuschaltung RK': ('Klimakälte', 'hoho'), \n",
    "    'Lastzuschaltung WW normiert Tag': ('Warmwasserbereitstellung', 'hoho')\n",
    "}\n",
    "\n",
    "# Change column names\n",
    "availability_ind_pos = map_column_names(availability_ind_pos, availability_cats_ind_pos)\n",
    "availability_ind_neg = map_column_names(availability_ind_neg, availability_cats_ind_neg)\n",
    "availability_tcs_pos = map_column_names(availability_tcs_pos, availability_cats_tcs_pos)\n",
    "availability_tcs_neg = map_column_names(availability_tcs_neg, availability_cats_tcs_neg)\n",
    "availability_hoho_pos = map_column_names(availability_hoho_pos, availability_cats_hoho_pos)\n",
    "availability_hoho_neg = map_column_names(availability_hoho_neg, availability_cats_hoho_neg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Introduce shortcut for availability time series\n",
    "availabilities_dict = {\n",
    "    (\"ind\", \"pos\"): availability_ind_pos,\n",
    "    (\"ind\", \"neg\"): availability_ind_neg,\n",
    "    (\"tcs\", \"pos\"): availability_tcs_pos,\n",
    "    (\"tcs\", \"neg\"): availability_tcs_neg,\n",
    "    (\"hoho\", \"pos\"): availability_hoho_pos,\n",
    "    (\"hoho\", \"neg\"): availability_hoho_neg\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create availability time series for demand response categories for which availability info is missing\n",
    "\n",
    "In general, there are three options for filling up the missing availability information. All of these are used in the order they are mentioned:\n",
    "1. Assign the existing availability time series of another process / demand response category due to large similarities (assumed).\n",
    "2. Create a synthetic load profile by defining availability factors for hours, weekdays and months and combining them (as done in Gils 2015).\n",
    "3. Assign a constant availability profile for the entire year when load is assumed to be rather time-invariant.\n",
    "\n",
    "The folowing availability information is used:\n",
    "* Industry sector:\n",
    "    * Foundries (German: \"Gießereien\") will be assigned the value for copper and zinc.\n",
    "    * Calcium carbide will be assigned the value for electric furnace steel.\n",
    "    * For the remaining categories, no profiles are available. As a first proxy, a constant availability profile is assumed.\n",
    "* Tcs sector:\n",
    "    * Process cold is assigned the value for storage cooling.\n",
    "    * In the first place, for all categories not covered, a constant availability profile is assumed. _&rarr; Needs to be changed!_\n",
    "* Household sector:\n",
    "    * fride-freezer combinations as well as freezers will be assigned the same profile as fridges.\n",
    "    * For all remaining categories not covered, a constant availability profile is assumed. _&rarr; Needs to be changed!_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Use availability time series of existing categories"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "is_executing": false,
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# Boadcast values from existing demand response categories\n",
    "\n",
    "# Industry sector\n",
    "ind_mapping = {\n",
    "    ('Calciumcarbid-Herstellung (Lichtbogenofen)', 'ind'): \n",
    "    ('Elektrostahlherstellung (Lichtbogenofen)', 'ind'),\n",
    "    ('Gießereien (Induktionsofen)', 'ind'): \n",
    "    ('Kupfer- und Zinkherstellung (Elektrolyse)', 'ind')\n",
    "}\n",
    "\n",
    "for k, v in ind_mapping.items():\n",
    "    availability_ind_pos[k] = availability_ind_pos[v]\n",
    "    availability_ind_neg[k] = availability_ind_neg[v]\n",
    "\n",
    "# tcs\n",
    "tcs_mapping = {\n",
    "    ('Prozesskälte', 'tcs'): ('Kühlhäuser', 'tcs')\n",
    "}\n",
    "\n",
    "for k, v in tcs_mapping.items():\n",
    "    availability_tcs_pos[k] = availability_tcs_pos[v]\n",
    "    availability_tcs_neg[k] = availability_tcs_neg[v]\n",
    "\n",
    "# households\n",
    "hoho_mapping = {\n",
    "    ('Kühl- und Gefrierkombinationen', 'hoho'): ('Kühlschränke', 'hoho'),\n",
    "    ('Gefrierschränke und -truhen', 'hoho'): ('Kühlschränke', 'hoho')\n",
    "}\n",
    "\n",
    "for k, v in hoho_mapping.items():\n",
    "    availability_hoho_pos[k] = availability_hoho_pos[v]\n",
    "    availability_hoho_neg[k] = availability_hoho_neg[v]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Determine for which categories availability time series are still missing\n",
    "for key, value in availabilities_dict.items():\n",
    "    determine_missing_cols(\n",
    "        dr_categories_unjoined.droplevel(2), \n",
    "        value, \n",
    "        sector=key[0], \n",
    "        direction=key[1]\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define synthetic profiles for availability\n",
    "Define **synthetic profiles** for the remaining demand response categories missing availability information.\n",
    "* These comprise of a factor for the hourly, daily (weekday information) and monthly patterns.\n",
    "* The three are multiplied in order to obtain hourly availability information.\n",
    "* A very similar approach is used in Gils (2015); assumptions are taken from Gils (2015), pp. 188-190 or own assumptions taken.\n",
    "* Positive potential is derived from assumed load pattern directly.\n",
    "* Negative potential is derived assuming a feasible load increase up to 1.2 the maximum value (or a monthly maximum >= 0 for seasonal dependency) and normalizing the whole thing again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Hourly factors\n",
    "hours = range(0, 24)\n",
    "\n",
    "hourly_factors_constant = {\n",
    "    \"pos\": [1] * 24,\n",
    "    \"neg\": [1] * 24\n",
    "}\n",
    "hourly_factors_morning_evening_reduced = {\n",
    "    \"pos\": [0.8] * 6 + [1.0] * 12 + [0.8] * 6,\n",
    "    \"neg\": [1] * 6 + [0.5] * 12 + [1] * 6\n",
    "}\n",
    "hourly_factors_day_reduced = {\n",
    "    \"pos\": [\n",
    "        0.85, 0.9, 1.0, 1.0, 1.0, 0.85,\n",
    "        0.7, 0.50, 0.50, 0.50, 0.50, 0.55,\n",
    "        0.55, 0.60, 0.65, 0.70, 0.75, 0.70, \n",
    "        0.60, 0.60, 0.80, 0.95, 0.95, 0.95\n",
    "    ],\n",
    "    \"neg\": [\n",
    "        0.50, 0.43, 0.29, 0.29, 0.29, 0.50, \n",
    "        0.71, 1.0, 1.0, 1.0, 1.0, 0.93, \n",
    "        0.86, 0.79, 0.71, 0.64, 0.71, 0.86, \n",
    "        0.86, 0.58, 0.57, 0.36, 0.36, 0.36\n",
    "    ]\n",
    "}\n",
    "        \n",
    "hourly_factors_climate_cold_ind = {\n",
    "    \"pos\": [\n",
    "        0.1, 0.1, 0.1, 0.2, 0.2, 0.2,\n",
    "        0.3, 0.4, 0.5, 0.6, 0.8, 1.0,\n",
    "        1.0, 1.0, 1.0, 1.0, 0.8, 0.7,\n",
    "        0.6, 0.5, 0.3, 0.2, 0.1, 0.1\n",
    "    ],\n",
    "    \"neg\": [\n",
    "        1.0, 1.0, 1.00, 0.91, 0.91, 0.91,\n",
    "        0.82, 0.73, 0.64, 0.55, 0.36, 0.18,\n",
    "        0.18, 0.18, 0.18, 0.18, 0.36, 0.45, \n",
    "        0.55, 0.64, 0.82, 0.91, 1.00, 1.00, \n",
    "    ]\n",
    "}\n",
    "\n",
    "hourly_factors = {\n",
    "    # constant\n",
    "    ('Prozesskälte', 'ind'): hourly_factors_constant,\n",
    "    ('Prozesswärme', 'ind'): hourly_factors_constant, \n",
    "    ('Luftzerlegung', 'ind'): hourly_factors_constant, \n",
    "    ('Prozesskälte', 'tcs'): hourly_factors_constant,\n",
    "    # morning & evening reduced\n",
    "    ('Druckluftanwendungen', 'ind'): hourly_factors_morning_evening_reduced,\n",
    "    ('Belüftung', 'ind'): hourly_factors_morning_evening_reduced,\n",
    "    ('Zerkleinerer', 'tcs'): hourly_factors_morning_evening_reduced,\n",
    "    ('Prozesswärme', 'tcs'): hourly_factors_morning_evening_reduced,\n",
    "    ('Heizungsumwälzpumpen', 'hoho'): hourly_factors_morning_evening_reduced,\n",
    "    # day reduced\n",
    "    ('Kühlung (Lebensmittelindustrie)', 'ind'): hourly_factors_day_reduced,\n",
    "    # special cases\n",
    "    ('Klimakälte', 'ind'): hourly_factors_climate_cold_ind,\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Weekly factors\n",
    "days = range(0, 7)\n",
    "\n",
    "weekly_factors_constant = {\n",
    "    \"pos\": [1] * 7,\n",
    "    \"neg\": [1] * 7\n",
    "}\n",
    "weekly_factors_weekend_slightly_reduced = {\n",
    "    \"pos\": [1] * 5 + [0.95, 0.9],\n",
    "    \"neg\": [0.67] * 5 + [0.83, 1]\n",
    "}\n",
    "weekly_factors_weekend_medium_reduced = {\n",
    "    \"pos\": [1] * 5 + [0.6, 0.5],\n",
    "    \"neg\": [0.29] * 5 + [0.86, 1] \n",
    "}\n",
    "weekly_factors_weekend_reduced = {\n",
    "    \"pos\": [1] * 5 + [0.5, 0.05],\n",
    "    \"neg\": [0.17] * 5 + [0.61, 1]\n",
    "}\n",
    "\n",
    "weekly_factors = {\n",
    "    # constant\n",
    "    ('Prozesskälte', 'ind'): weekly_factors_constant,\n",
    "    ('Prozesswärme', 'ind'): weekly_factors_constant,\n",
    "    ('Luftzerlegung', 'ind'): weekly_factors_constant,\n",
    "    ('Heizungsumwälzpumpen', 'hoho'): weekly_factors_constant,\n",
    "    # weekend sligthly reduced\n",
    "    ('Klimakälte', 'ind'): weekly_factors_weekend_slightly_reduced,\n",
    "    ('Kühlung (Lebensmittelindustrie)', 'ind'): weekly_factors_weekend_slightly_reduced,\n",
    "    ('Druckluftanwendungen', 'ind'): weekly_factors_weekend_slightly_reduced,\n",
    "    ('Zerkleinerer', 'tcs'): weekly_factors_weekend_slightly_reduced,\n",
    "    # weekend medium reduced\n",
    "    ('Belüftung', 'ind'): weekly_factors_weekend_medium_reduced,\n",
    "    # weekend reduced\n",
    "    ('Prozesskälte', 'tcs'): weekly_factors_weekend_reduced,\n",
    "    ('Prozesswärme', 'tcs'): weekly_factors_weekend_reduced,\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Monthly factors\n",
    "months = range(1, 13)\n",
    "\n",
    "monthly_factors_constant = {\n",
    "    \"pos\": [1] * 12,\n",
    "    \"neg\": [1] * 12\n",
    "}\n",
    "monthly_factors_reduced_in_winter = {\n",
    "    \"pos\": [0.9] * 2 + [0.95] + [1] * 6 + [0.95] + [0.9] * 2,\n",
    "    \"neg\": [0.5] * 2 + [0.75] + [1] * 6 + [0.75] + [0.5] * 2\n",
    "}\n",
    "monthly_factors_warm_seasons = {\n",
    "    \"pos\": [0] * 4 + [0.3, 0.7] + [1] * 2 + [0.6, 0.1] + [0] * 2,\n",
    "    \"neg\": [0] * 4 + [1.0, 0.5] + [0.5] * 2 + [0.5, 1.0] + [0] * 2\n",
    "}\n",
    "monthly_factors_heating_seasons = {\n",
    "    \"pos\": [1, 0.8, 0.5, 0.1] + 5 * [0] + [0.2, 0.5, 1],\n",
    "    \"neg\": [0.5, 0.5, 0.75, 1.0] + 5 * [0] + [1.0, 0.75, 0.5]\n",
    "}\n",
    "\n",
    "monthly_factors = {\n",
    "    # constant\n",
    "    ('Prozesskälte', 'ind'): monthly_factors_constant,\n",
    "    ('Prozesswärme', 'ind'): monthly_factors_constant,\n",
    "    ('Luftzerlegung', 'ind'): monthly_factors_constant,\n",
    "    ('Druckluftanwendungen', 'ind'): monthly_factors_constant,\n",
    "    ('Belüftung', 'ind'): monthly_factors_constant,\n",
    "    ('Zerkleinerer', 'tcs'): monthly_factors_constant,\n",
    "    ('Prozesswärme', 'tcs'): monthly_factors_constant,\n",
    "    # reduced in winter\n",
    "    ('Kühlung (Lebensmittelindustrie)', 'ind'): monthly_factors_reduced_in_winter,\n",
    "    ('Prozesskälte', 'tcs'): monthly_factors_reduced_in_winter,\n",
    "    # warm seasons\n",
    "    ('Klimakälte', 'ind'): monthly_factors_warm_seasons,\n",
    "    # heating seasons\n",
    "    ('Heizungsumwälzpumpen', 'hoho'): monthly_factors_heating_seasons\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "factors = create_synthetic_profile_factors(\n",
    "    {\n",
    "        hours: (\"hours\", hourly_factors),\n",
    "        days: (\"days\", weekly_factors),\n",
    "        months: (\"months\", monthly_factors)\n",
    "    }\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Determine process categories per sector that need to be assigned synthetic profiles\n",
    "synthetic_cols = dict()\n",
    "for sector in sectors:\n",
    "    synthetic_cols[sector] = [col for col in hourly_factors.keys() if col[1] == sector]\n",
    "\n",
    "# Assign synthetic profiles per sector and direction\n",
    "for key, value in availabilities_dict.items():\n",
    "    sector = key[0]\n",
    "    direction = key[1]\n",
    "    \n",
    "    # No missing availability for hoho pos -> skip!\n",
    "    if not key == (\"hoho\", \"pos\"):\n",
    "        assign_availability_remaining(\n",
    "            parameters_for_clustering_unjoined, \n",
    "            value,\n",
    "            synthetic_cols=synthetic_cols[sector],\n",
    "            factors=factors,\n",
    "            sector=sector,\n",
    "            direction=direction\n",
    "        )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create load profiles per cluster\n",
    "Create load profiles for the demand response categories.\n",
    "In principle, two approaches could be used here:\n",
    "1. Most simple proxy: Use availability time series in positive direction. &rarr; drawback: very rough estimate\n",
    "2. More advanced: Use profiles per WZ as an intermediate product of the the demand regio disaggregator tool.\n",
    "\n",
    "Profiles are scaled with installed capacity in both cases. Simultaneity factors are introduced to account for the fact that not all demands happen simultaneously."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Preparation for load profile assessment\n",
    "* Provide availability time series with cluster information\n",
    "* Initialize empty dictionaries to store the time series\n",
    "* Determine which max. simultaneity factors to use &rarr; assumption is needed since there is no complete information on energy consumption or full load hours for every appliance\n",
    "* Copy availability time series in order to be able to use it for both approaches\n",
    "\n",
    "_NOTE: One could derive simultaneity from a synthetic \"Gleichzeitigkeitsfunktion\" when full load hours were given ..._"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Assign cluster information\n",
    "for key in availabilities_dict.keys():\n",
    "    availabilities_dict[key].loc[\"cluster\"] = parameters_for_clustering_unjoined.loc[\n",
    "        parameters_for_clustering_unjoined[\"year\"] == \"SQ\", \"cluster\"\n",
    "    ]\n",
    "    if key in [(\"tcs\", \"pos\"), (\"tcs\", \"neg\")]:\n",
    "        # Assign missing cluster information for heat pumps (solely attributed to households)\n",
    "        availabilities_dict[key][(\"Wärmepumpen\", \"tcs\")].loc[\"cluster\"] = (\n",
    "            parameters_for_clustering_unjoined.loc[\n",
    "                parameters_for_clustering_unjoined[\"year\"] == \"SQ\"\n",
    "            ].at[(\"Wärmepumpen\", \"hoho\"), \"cluster\"]\n",
    "        )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "timeseries_dict = dict()\n",
    "installed_cap_dict = dict()\n",
    "\n",
    "# determine maximum simultaneity factors per sector\n",
    "max_simultaneity = {\n",
    "    \"ind\": 0.9,\n",
    "    \"tcs\": 0.4,\n",
    "    \"hoho\": 0.1\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a copy of availability time series\n",
    "load_profile_ind = availability_ind_pos.copy()\n",
    "load_profile_tcs = availability_tcs_pos.copy()\n",
    "load_profile_hoho = availability_hoho_pos.copy()\n",
    "\n",
    "# Create a dict shortcut to access\n",
    "load_profiles_dict = {\n",
    "    \"ind\": load_profile_ind,\n",
    "    \"tcs\": load_profile_tcs,\n",
    "    \"hoho\": load_profile_hoho\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Approach 1: Use data from the demand regio disaggregator\n",
    "* Read in normalized profiles as output from the disaggregator tool\n",
    "* Map the demand response profiles with the \"Wirtschaftszweige\" (WZ) used in the disaggregator tool\n",
    "    * The mapping has to be done / updated manually by selecting the appropriate / closest WZ from the DESTATIS categorization.\n",
    "    * Households are assigned a value of 0 (corresponding to the usage of the standard load profile H0).\n",
    "    * Cross-cutting technologies are assigned a value of -1.\n",
    "* For the WZs which can be directly assigned, the respective demand pattern is used\n",
    "* For cross-cutting technologies no data exists, hence, the availability time series is used as a backup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not use_ava_ts_for_profiles:\n",
    "    profiles_ind_normalized = pd.read_csv(\n",
    "        path_folder_in+\"profiles_industry_normalized.csv\", \n",
    "        index_col=0, \n",
    "        parse_dates=True\n",
    "    )\n",
    "    profiles_tcs_normalized = pd.read_csv(\n",
    "        path_folder_in+\"profiles_cts_normalized.csv\", \n",
    "        index_col=0, \n",
    "        parse_dates=True\n",
    "    )\n",
    "    profiles_hoho_normalized = pd.read_csv(\n",
    "        path_folder_in+\"profiles_households_normalized.csv\",\n",
    "        index_col=0,\n",
    "        parse_dates=True\n",
    "    )\n",
    "\n",
    "    WZ_mapping = pd.read_csv(\n",
    "        path_folder_in+\"remaining_categories_WZ_mapping.csv\", \n",
    "        index_col=[0, 1],\n",
    "        sep=\";\"\n",
    "    )\n",
    "\n",
    "    profiles_ind_normalized.columns = profiles_ind_normalized.columns.astype(int)\n",
    "    profiles_tcs_normalized.columns = profiles_tcs_normalized.columns.astype(int)\n",
    "    profiles_hoho_normalized.columns = profiles_hoho_normalized.columns.astype(int)\n",
    "    \n",
    "    # Prepare WZ mapping\n",
    "    WZ_mapping.index.names = [\"Prozesskategorie\", \"help_sector\"]\n",
    "    WZ_mapping[\"sector\"] = np.where(\n",
    "        WZ_mapping.index.get_level_values(1) == \"tcs+hoho\", \"hoho\",\n",
    "        WZ_mapping.index.get_level_values(1)\n",
    "    )\n",
    "    \n",
    "    WZ_mapping = WZ_mapping.reset_index(level=1)\n",
    "    WZ_mapping = WZ_mapping.set_index(\"sector\", append=True)\n",
    "\n",
    "    WZs_used = WZ_mapping[\"WZ\"].unique()\n",
    "    \n",
    "    display(WZ_mapping.head(10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not use_ava_ts_for_profiles:\n",
    "    normalized_profiles = {\n",
    "        \"ind\": profiles_ind_normalized,\n",
    "        \"tcs\": profiles_tcs_normalized,\n",
    "        \"hoho\": profiles_hoho_normalized\n",
    "    }\n",
    "\n",
    "    load_timeseries_dict = dict()\n",
    "    processes_used_dict = {\"ind\": [], \"tcs\": [], \"hoho\": []}\n",
    "    all_processes_used = []\n",
    "\n",
    "    for sector in normalized_profiles.keys():\n",
    "\n",
    "        load_timeseries_dict[sector] = pd.DataFrame(\n",
    "            index=normalized_profiles[sector].index,\n",
    "            columns=parameters_for_clustering_unjoined.index[\n",
    "                (parameters_for_clustering_unjoined[\"year\"] == \"SQ\")\n",
    "                & (parameters_for_clustering_unjoined.index.get_level_values(1) == sector)\n",
    "            ]\n",
    "        )\n",
    "\n",
    "\n",
    "    for WZ in WZs_used:\n",
    "        # Get processes for which WZ can be used\n",
    "        ix = [process for process in WZ_mapping[WZ_mapping[\"WZ\"] == WZ].index]\n",
    "        \n",
    "        for sector, profile in normalized_profiles.items():\n",
    "            if WZ in normalized_profiles[sector].columns:\n",
    "                # Assign values of respective WZ\n",
    "                reshaped_vals = np.reshape(\n",
    "                    profile[WZ].values,\n",
    "                    [len(profile[WZ].values),1]\n",
    "                )\n",
    "                reshaped_vals = np.repeat(reshaped_vals, len(ix), axis=1)\n",
    "                \n",
    "                load_timeseries_dict[sector][ix] = reshaped_vals\n",
    "                \n",
    "                processes_used_dict[sector].extend(ix)\n",
    "                all_processes_used.extend(ix)\n",
    "\n",
    "    remaining_processes = [\n",
    "        ix for ix in parameters_for_clustering_unjoined[\n",
    "            parameters_for_clustering_unjoined[\"year\"] == \"SQ\"\n",
    "        ].index\n",
    "        if ix not in all_processes_used\n",
    "    ]\n",
    "    \n",
    "    for sector, profile in load_profiles_dict.items():\n",
    "        # Replace the default load time series with load profile values and (re-)add cluster info\n",
    "        profile[processes_used_dict[sector]] = load_timeseries_dict[\n",
    "            sector\n",
    "        ][processes_used_dict[sector]]\n",
    "        profile.loc[\"cluster\"] = availabilities_dict[(sector, \"pos\")].loc[\"cluster\"]\n",
    "    \n",
    "    # Print logging on which proxy to use\n",
    "    dr_proxy_string =  (\n",
    "        \"Using branch profile from demand regio disaggregator \"\n",
    "        \"as a load profile proxy for the following categories:\"\n",
    "    )\n",
    "    print(dr_proxy_string)\n",
    "    print(\"-\" * len(dr_proxy_string))\n",
    "    for process in all_processes_used:\n",
    "        print(process)\n",
    "    \n",
    "    print()\n",
    "\n",
    "    ava_proxy_string = (\n",
    "        \"Using positive availability time series as a load profile proxy for the following categories:\"\n",
    "    )\n",
    "    print(ava_proxy_string)\n",
    "    print(\"-\" * len(ava_proxy_string))\n",
    "    for process in remaining_processes:\n",
    "        print(process)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Approach 2: Use positive availability time series as simplest proxy\n",
    "> _NOTE: If approach 1 is chosen, this will be run as well, but the data used will differ!_\n",
    "* Multiply availability time series with installed capacity and maximum simultaneity factor for the respective sector\n",
    "* Include cluster information for grouping.\n",
    "* Combine to an overall load profile time series dataset which will be grouped within the clusters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for year in years:\n",
    "    for col in quantile_cols:\n",
    "        \n",
    "        to_concat = []\n",
    "\n",
    "        for sector in load_profiles_dict.keys():\n",
    "            installed_cap_dict[(sector, year, col)] = parameters_for_clustering_unjoined[\n",
    "                (parameters_for_clustering_unjoined[\"year\"] == year) \n",
    "                & (parameters_for_clustering_unjoined.index.get_level_values(1) == sector)\n",
    "            ][col+\"_installed_cap\"]\n",
    "            \n",
    "            timeseries_dict[(sector, year, col)] = pd.DataFrame(\n",
    "                index=load_profiles_dict[sector].index, \n",
    "                columns=load_profiles_dict[sector].columns\n",
    "            )\n",
    "            \n",
    "            timeseries_dict[(sector, year, col)] = load_profiles_dict[sector].iloc[:-1].mul(\n",
    "                installed_cap_dict[(sector, year, col)]).mul(max_simultaneity[sector])\n",
    "\n",
    "            timeseries_dict[(sector, year, col)].loc[\"cluster\"] = load_profiles_dict[sector].loc[\"cluster\"]\n",
    "            \n",
    "            to_concat.append(timeseries_dict[(sector, year, col)])\n",
    "        \n",
    "        timeseries_dict[(\"all_sectors\", year, col)] = pd.concat(to_concat, axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualize load profiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(15, 5))\n",
    "timeseries_dict[(\"ind\", \"SQ\", \"50%\")].iloc[:672].plot(ax=ax)\n",
    "plt.ylim([0,2500])\n",
    "plt.legend(bbox_to_anchor=(1.05, 1))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(15, 5))\n",
    "timeseries_dict[(\"tcs\", \"SQ\", \"50%\")].iloc[:672].plot(ax=ax)\n",
    "plt.ylim([0,1500])\n",
    "plt.legend(bbox_to_anchor=(1.05, 1))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(15, 5))\n",
    "timeseries_dict[(\"hoho\", \"SQ\", \"50%\")].iloc[:672].plot(ax=ax)\n",
    "plt.ylim([0,7000])\n",
    "plt.legend(bbox_to_anchor=(1.05, 1))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create availability time series per cluster\n",
    "* Calculate weighted averages within clusters:\n",
    "    * Iterate over all years and clusters\n",
    "    * Determine simultaneous maximum positive and negative potentials taking into account eligibility for shifting, shedding or both\n",
    "    * Assign simultaneous values as summed potentials\n",
    "    * Rescale to a value of 1\n",
    "* Calculate summed load profiles as well\n",
    "* Put all results into dict structures\n",
    "* Save the results, i.e. availability time series per cluster, to a csv file."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Interpretation of availability factors per cluster**:\n",
    "* Maximum availability per cluster does not necessarily have to be 1.\n",
    "* A value of 0.88 implies that at maximum 88% of the overall cluster capacity are available at the same time. &rarr; Similar interpretation than simultaneity factor.\n",
    "* Nevertheless, a rescaling is done here for the sake of easier interpretation. &rarr; I.e. maximum values will be 1.0 and cluster capacity is adapted (scaled_down) accordingly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sector_dict = {}\n",
    "\n",
    "for year in years:\n",
    "    \n",
    "    # Sector dict contains: original availabilites (pos & neg) and parameter datas\n",
    "    sector_dict[year] = {\n",
    "        \"ind\": (\n",
    "            availability_ind_pos, \n",
    "            availability_ind_neg, \n",
    "            parameters_for_clustering_unjoined.loc[\n",
    "                (parameters_for_clustering_unjoined[\"year\"] == year)\n",
    "                & (parameters_for_clustering_unjoined.index.get_level_values(1) == \"ind\")\n",
    "            ]\n",
    "        ),\n",
    "        \"tcs\": (\n",
    "            availability_tcs_pos, \n",
    "            availability_tcs_neg, \n",
    "            parameters_for_clustering_unjoined.loc[\n",
    "                (parameters_for_clustering_unjoined[\"year\"] == year)\n",
    "                & (parameters_for_clustering_unjoined.index.get_level_values(1) == \"tcs\")\n",
    "            ]\n",
    "        ),\n",
    "        \"hoho\": (\n",
    "            availability_hoho_pos, \n",
    "            availability_hoho_neg, \n",
    "            parameters_for_clustering_unjoined.loc[\n",
    "                (parameters_for_clustering_unjoined[\"year\"] == year)\n",
    "                & (parameters_for_clustering_unjoined.index.get_level_values(1) == \"hoho\")\n",
    "            ]\n",
    "        )\n",
    "    }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract info on combined clusters for tcs and hoho separately in order not to overwrite it\n",
    "tcs_hoho_clusters = parameters_for_clustering.loc[\n",
    "    (parameters_for_clustering[\"year\"] == \"SQ\")\n",
    "    & (parameters_for_clustering.index.get_level_values(1) == \"tcs+hoho\"),\n",
    "    \"cluster\"\n",
    "].unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Renaming for nicer plots\n",
    "sector_map = {\n",
    "    \"hoho\": \"Haushalte\",\n",
    "    \"tcs\": \"GHD\",\n",
    "    \"ind\": \"Industrie\",\n",
    "    \"tcs+hoho\": \"Haushalte und GHD\",\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pot_cols = [\"potential_pos_overall\", \"potential_neg_overall\", \"installed_cap\"]\n",
    "cols_of_interest = [i+\"_\"+j for i in quantile_cols for j in pot_cols]\n",
    "\n",
    "# Use dicts and DataFrames to store overall cluster potential and availability time series per cluster\n",
    "cluster_overall_pot_dict = {}\n",
    "cluster_overall_ts_dict = {}\n",
    "availability_clusters = pd.DataFrame(index=availability_ind_pos.index).drop(\"cluster\")\n",
    "\n",
    "# Create availability time series within clusters by calculating weighted averages\n",
    "for year in years:\n",
    "    \n",
    "    cluster_overall_pot_df = pd.DataFrame(columns=cols_of_interest)\n",
    "    load_timeseries = pd.DataFrame(index=timeseries_dict[(\"ind\", \"SQ\",\"50%\")].iloc[:-1].index)\n",
    "    \n",
    "    for cluster_number in np.sort(parameters_for_clustering[\"cluster\"].unique()):\n",
    "        \n",
    "        to_concat = []\n",
    "        cluster_data = {}\n",
    "        \n",
    "        for col in quantile_cols:\n",
    "            sector = parameters_for_clustering_unjoined.loc[\n",
    "                (parameters_for_clustering_unjoined[\"year\"] == year)\n",
    "                & (parameters_for_clustering_unjoined[\"cluster\"] == cluster_number)\n",
    "            ].index.get_level_values(1)[0]\n",
    "            \n",
    "            # introduce shortcuts for readability\n",
    "            # original availability series and extract from parameter data set\n",
    "            org_ava_pos = sector_dict[year][sector][0]\n",
    "            org_ava_neg = sector_dict[year][sector][1]\n",
    "            potentials = sector_dict[year][sector][2]\n",
    "            \n",
    "            # Calculate a weighted average for positive potentials \n",
    "            # (weights: maximum overall potential information)\n",
    "            ava_pos = org_ava_pos.loc[\n",
    "                :, org_ava_pos.loc[\"cluster\"] == cluster_number\n",
    "            ].drop(\"cluster\")\n",
    "            \n",
    "            pot_pos = potentials.loc[\n",
    "                potentials[\"cluster\"] == cluster_number\n",
    "            ]\n",
    "            kind = parameters_for_clustering_unjoined.loc[\n",
    "                (parameters_for_clustering_unjoined[\"year\"] == year)\n",
    "                & (parameters_for_clustering_unjoined[\"cluster\"] == cluster_number)\n",
    "            ][\"kind\"].unique()[0]\n",
    "            \n",
    "            # Choose positive potential columns dependent on eligibility:\n",
    "            # shifting only: regular positive potential\n",
    "            # shedding only: positive shedding potential\n",
    "            # shifting and shedding: element-wise maximum of both (since they are competing)\n",
    "            pot_pos_cols = [\"_potential_pos_overall\"]\n",
    "            if kind == \"shed_only\":\n",
    "                pot_pos_cols = [\"_potential_pos_overall_shed\"]\n",
    "            if kind == \"shift_shed\":\n",
    "                pot_pos_cols.append(\"_potential_pos_overall_shed\")\n",
    "            \n",
    "            pot_pos = pot_pos[[\n",
    "                col + pot_pos_col \n",
    "                for pot_pos_col in pot_pos_cols\n",
    "            ]].max(axis=1)\n",
    "            \n",
    "            cluster_overall_pot_pos = pot_pos.sum()\n",
    "            \n",
    "            if cluster_overall_pot_pos != 0:\n",
    "                ava_pos[col+\"_pos_cluster_\"+str(cluster_number)] = ava_pos.mul(\n",
    "                    pot_pos\n",
    "                ).div(cluster_overall_pot_pos).sum(axis=1)\n",
    "\n",
    "                # Do rescaling: Adjust (reduce) potential information if necessary \n",
    "                max_pot_pos = (\n",
    "                    ava_pos[col+\"_pos_cluster_\"+str(cluster_number)].max() \n",
    "                    * cluster_overall_pot_pos\n",
    "                )\n",
    "                # Scale max value of availability time series to 1   \n",
    "                ava_pos[col+\"_pos_cluster_\"+str(cluster_number)] = (\n",
    "                    ava_pos[col+\"_pos_cluster_\"+str(cluster_number)].div(\n",
    "                        ava_pos[col+\"_pos_cluster_\"+str(cluster_number)].max()\n",
    "                    )\n",
    "                )\n",
    "\n",
    "            # Calculate a weighted average for negative potentials \n",
    "            # (weights: maximum overall potential information)\n",
    "            ava_neg = org_ava_neg.loc[\n",
    "                :, org_ava_neg.loc[\"cluster\"] == cluster_number\n",
    "            ].drop(\"cluster\")\n",
    "            \n",
    "            # Choose negative potential dependent on eligibility:\n",
    "            # shedding only: No negative potential existing\n",
    "            if kind == \"shed_only\":\n",
    "                pot_neg = 0\n",
    "                cluster_overall_pot_neg = 0\n",
    "            else:\n",
    "                pot_neg = potentials.loc[\n",
    "                    potentials[\"cluster\"] == cluster_number, \n",
    "                    col+\"_potential_neg_overall\"\n",
    "                ]\n",
    "                cluster_overall_pot_neg = pot_neg.sum()\n",
    "            \n",
    "            if cluster_overall_pot_neg != 0:\n",
    "                ava_neg[col+\"_neg_cluster_\"+str(cluster_number)] = ava_neg.mul(\n",
    "                    pot_neg\n",
    "                ).div(cluster_overall_pot_neg).sum(axis=1)\n",
    "\n",
    "                # Do rescaling: Adjust (reduce) potential information if necessary \n",
    "                max_pot_neg = (\n",
    "                    ava_neg[col+\"_neg_cluster_\"+str(cluster_number)].max() \n",
    "                    * cluster_overall_pot_neg\n",
    "                )\n",
    "                # Scale max value of availability time series to 1   \n",
    "                ava_neg[col+\"_neg_cluster_\"+str(cluster_number)] = (\n",
    "                    ava_neg[col+\"_neg_cluster_\"+str(cluster_number)].div(\n",
    "                        ava_neg[col+\"_neg_cluster_\"+str(cluster_number)].max()\n",
    "                    )\n",
    "                )\n",
    "            else:\n",
    "                ava_neg[col+\"_neg_cluster_\"+str(cluster_number)] = 0\n",
    "                max_pot_neg = 0\n",
    "            \n",
    "            # Use proper cluster number label (int instead of float)\n",
    "            if cluster_number not in tcs_hoho_clusters:\n",
    "                cluster_label = sector+\"_cluster-\"+str(int(cluster_number))\n",
    "            else:\n",
    "                cluster_label = \"tcs+hoho_cluster-\"+str(int(cluster_number))\n",
    "\n",
    "            # show exemplarily availability time series for clusters\n",
    "            if year == \"SQ\":\n",
    "                availability_clusters[col+\"_\"+cluster_label+\"_pos\"] = (\n",
    "                    ava_pos[col+\"_pos_cluster_\"+str(cluster_number)]\n",
    "                )\n",
    "                availability_clusters[col+\"_\"+cluster_label+\"_neg\"] = (\n",
    "                    ava_neg[col+\"_neg_cluster_\"+str(cluster_number)]\n",
    "                )\n",
    "                availability_clusters = availability_clusters.round(4)\n",
    "            \n",
    "                if col == \"50%\":\n",
    "                    colors.at[\"Mittelwert Cluster\", \"Farbe (matplotlib strings)\"] = \"darkblue\"\n",
    "                    ava_pos_plot = ava_pos.copy()\n",
    "                    if cluster_number == 101:\n",
    "                        ava_pos_plot = ava_pos_plot.rename(columns={col[1]: \"tcs+hoho\" for col in ava_pos_plot.columns if not \"cluster\" in col[0]})\n",
    "                    ava_pos_plot = ava_pos_plot.rename(columns={col[1]: sector_map[col[1]] for col in ava_pos_plot.columns if not \"cluster\" in col[0]})    \n",
    "                    ava_pos_plot.columns = ava_pos_plot.columns.get_level_values(0)\n",
    "                    ava_pos_plot = ava_pos_plot.rename(columns={col: \"Mittelwert Cluster\" for col in ava_pos_plot.columns if \"cluster\" in col})\n",
    "                    ava_pos_plot.index = pd.date_range(start=\"2017-01-01 00:00:00\", periods=len(ava_pos), freq=\"H\")\n",
    "                    \n",
    "                    # Show a sample week for the clusters \n",
    "                    # ... for status quo, positive potentials and median values only\n",
    "                    fig, ax = plt.subplots(figsize=(12, 6))\n",
    "                    _ = ava_pos_plot.iloc[169:362,:-1].plot(ax=ax, color=colors[\"Farbe (matplotlib strings)\"].to_dict(), linewidth=2)\n",
    "                    _ = ava_pos_plot.iloc[169:362,-1:].plot(ax=ax, linewidth=5, color=colors[\"Farbe (matplotlib strings)\"].to_dict())\n",
    "                    #_ = plt.title(\"Cluster \"+str(cluster_number))\n",
    "                    _ = plt.ylim([-0.05, 1.05])\n",
    "                    _ = plt.legend(loc=\"upper center\", bbox_to_anchor=(0.5, -0.28), fancybox=True, ncol=4)\n",
    "                    _ = plt.xlabel(\"Zeit\", labelpad=10)\n",
    "                    _ = plt.ylabel(\"Normiertes Lastreduktionspotenzial\", labelpad=10)\n",
    "                    _ = plt.tight_layout()\n",
    "                    current_values = plt.gca().get_yticks()\n",
    "                    _ = plt.gca().set_yticklabels(['{:,.2f}'.format(x).replace(\".\", \",\") for x in current_values])\n",
    "                    plt.show()\n",
    "                    fig.savefig(path_folder_plots+\"cluster_\"+str(cluster_number)+\".png\", dpi=300, bbox_inches=\"tight\")\n",
    "                    plt.close()\n",
    "            \n",
    "            # Save potential outputs\n",
    "            cluster_data[col+\"_potential_pos_overall\"] = max_pot_pos\n",
    "            cluster_data[col+\"_potential_neg_overall\"] = max_pot_neg\n",
    "            \n",
    "            # Assign load profiles and store them in a dict\n",
    "            load_timeseries[col+\"_\"+cluster_label] = (\n",
    "                timeseries_dict[(\"all_sectors\", year, col)].loc[\n",
    "                    :, timeseries_dict[(\"all_sectors\", year, col)].loc[\"cluster\"] \n",
    "                    == cluster_number\n",
    "                ].sum(axis=1).drop(\"cluster\")\n",
    "            )\n",
    "            \n",
    "            max_cap = load_timeseries[col+\"_\"+cluster_label].max()\n",
    "            \n",
    "            # Save installed cluster capacity\n",
    "            cluster_data[col+\"_installed_cap\"] = max_cap\n",
    "            \n",
    "        # Combine potential outputs and store them in dict of DataFrames\n",
    "        cluster_overall_pot_df.loc[\n",
    "            cluster_label, [\n",
    "                i+\"_\"+j for i in quantile_cols \n",
    "                for j in pot_cols\n",
    "            ]] = cluster_data\n",
    "\n",
    "        cluster_overall_pot_df.loc[\n",
    "            cluster_label, \"kind\"\n",
    "        ] = kind\n",
    "        \n",
    "        cluster_overall_pot_dict[year] = cluster_overall_pot_df\n",
    "        cluster_overall_ts_dict[year] = load_timeseries\n",
    "\n",
    "display(availability_clusters.tail())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Introduce some fixes:\n",
    "* Normalize load profiles again in order to be able to use them combined with maximum capacity demand\n",
    "* Rename `installed_cap` to `max_cap` in order to prevent misinterpretation\n",
    "* Limit positive potential to the value of the respective maximum simultaneously dispatched capacity, i.e. `max_cap`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Normalize load profiles (for usage with pommes / oemof-solph)\n",
    "# Rename column to 'max_cap' in order to be able to distinct it from installed capacity\n",
    "for year in years:\n",
    "    cluster_overall_ts_dict[year] = cluster_overall_ts_dict[year].div(\n",
    "        cluster_overall_ts_dict[year].max()\n",
    "    )\n",
    "    for col in quantile_cols:\n",
    "        cluster_overall_pot_dict[year].rename(\n",
    "            {col + \"_installed_cap\": col + \"_max_cap\" for col in quantile_cols},\n",
    "            axis=1,\n",
    "            inplace=True\n",
    "        )\n",
    "    \n",
    "    for col in quantile_cols:\n",
    "        cluster_overall_pot_dict[year][col+\"_potential_pos_overall\"] = (\n",
    "            cluster_overall_pot_dict[year][\n",
    "                [col+\"_potential_pos_overall\", col+\"_max_cap\"]\n",
    "            ].min(axis=1)\n",
    "        )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cluster_overall_pot_dict[\"SQ\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cluster_overall_pot_dict[\"2050\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cluster_overall_ts_dict[\"2050\"].max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cluster_labels = cluster_overall_pot_dict[\"SQ\"].index.values\n",
    "cluster_ts_by_cols = {}\n",
    "ava_pos_ts_by_cols = {}\n",
    "ava_neg_ts_by_cols = {}\n",
    "\n",
    "# Split timeseries into subsets\n",
    "for col in quantile_cols:\n",
    "    for year in years:\n",
    "        ava_cols = [i+\"_\"+j for i in [col] for j in cluster_labels]\n",
    "        ava_cols_pos = [i+\"_\"+j+\"_pos\" for i in [col] for j in cluster_labels]\n",
    "        ava_cols_neg = [i+\"_\"+j+\"_neg\" for i in [col] for j in cluster_labels]\n",
    "\n",
    "        cluster_ts_by_cols[col+\"_\"+year] = cluster_overall_ts_dict[year][ava_cols]\n",
    "        ava_pos_ts_by_cols[col+\"_\"+year] = availability_clusters[ava_cols_pos]\n",
    "        ava_neg_ts_by_cols[col+\"_\"+year] = availability_clusters[ava_cols_neg]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if write_outputs:\n",
    "    availability_clusters.to_csv(\n",
    "        path_folder_availability+filename_availability_out, \n",
    "        sep=\";\", \n",
    "        decimal=\",\"\n",
    "    )\n",
    "    write_multiple_sheets(\n",
    "        cluster_ts_by_cols, \n",
    "        path_folder_parameterization, \n",
    "        filename_load_profiles_out+\".xlsx\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Group the data within the clusters and write outputs\n",
    "* Determine, how the grouping will take place and which aggregation rules to apply for a parameter by specifying them in a dictionary\n",
    "* Perform the actual aggregation\n",
    "* Rename the clusters to a more human readable form and do some slight renamings\n",
    "* Add an assumption on unit lifetime and efficiency"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract parameters of interest for further model-based analyses\n",
    "params_of_interest = params_to_use + params_remaining\n",
    "params_of_interest = [param[1] for param in params_of_interest]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract aggregation rules\n",
    "grouping_cols = [\"cluster\"]\n",
    "agg_rules = {}\n",
    "\n",
    "for param, agg_rule in parameters_agg_dict.items():\n",
    "    if param not in params_of_interest:\n",
    "        continue\n",
    "    else:\n",
    "        for col in quantile_cols:\n",
    "            agg_rules[(col, param)] = agg_rule"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for col in quantile_cols:\n",
    "    agg_rules[col] = {\n",
    "        key[0] + \"_\" + key[1]: value for key, value in agg_rules.items()\n",
    "        if col == key[0]\n",
    "    }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Rename clusters to more intuitive names\n",
    "name_dict = {\n",
    "    \"tcs+hoho_cluster-101\": \"tcs+hoho_cluster_shift_only\",\n",
    "    \"hoho_cluster-201\": \"hoho_cluster_shift_shed\",\n",
    "    \"hoho_cluster-301\": \"hoho_cluster_shift_only\",\n",
    "    \"tcs_cluster-401\": \"tcs_cluster_shift_only\",\n",
    "    \"ind_cluster-601\": \"ind_cluster_shed_only\",\n",
    "    \"ind_cluster-701\": \"ind_cluster_shift_shed\",\n",
    "    \"ind_cluster-801\": \"ind_cluster_shift_only\"\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the parameters for the clusters\n",
    "overall_dict = {}\n",
    "pot_cols = [\n",
    "    \"potential_pos_overall\",\n",
    "    \"potential_neg_overall\",\n",
    "    \"max_cap\"\n",
    "]\n",
    "dur_cols = [\n",
    "    \"interference_duration_neg\",\n",
    "    \"interference_duration_pos\", \n",
    "    \"interference_duration_pos_shed\", \n",
    "    \"regeneration_duration\", \n",
    "    \"shifting_duration\"\n",
    "]\n",
    "cost_cols = [\n",
    "    \"fixed_costs\",\n",
    "    \"variable_costs\",\n",
    "    \"variable_costs_shed\",\n",
    "    \"specific_investments\"\n",
    "]\n",
    "other_cols = [\n",
    "    \"maximum_activations_year\",\n",
    "    \"installed_cap\", \"max_cap\",\n",
    "    \"potential_neg_overall\",\n",
    "    \"potential_pos_overall\"\n",
    "]\n",
    "\n",
    "parameters_for_clustering = parameters_for_clustering.set_index(\"kind\", append=True)\n",
    "\n",
    "for year in years:\n",
    "\n",
    "    params_year = parameters_for_clustering.loc[\n",
    "        parameters_for_clustering[\"year\"] == year\n",
    "    ].drop(columns=\"year\").astype(\"float64\").reset_index(level=[1, 2])\n",
    "        \n",
    "    for col in quantile_cols:\n",
    "        \n",
    "        overall_dict[col+\"_\"+year] = group_potential(\n",
    "            params_year, \n",
    "            grouping_cols=[\"cluster\", \"sector\", \"kind\"],\n",
    "            agg_rules=agg_rules[col],\n",
    "            # potential pos overall is used due to it giving the best data basis\n",
    "            weight_col=col+\"_potential_pos_overall\"\n",
    "        ).round(2)\n",
    "        \n",
    "        cols_potentials = [i+\"_\"+j for i in [col] for j in pot_cols]\n",
    "\n",
    "        # Update potential data with availability information from above\n",
    "        if adjust_potentials:\n",
    "            overall_dict[col+\"_\"+year][cols_potentials] = (\n",
    "                cluster_overall_pot_dict[year][cols_potentials].astype(float).round(2)\n",
    "            )\n",
    "\n",
    "        # Add country and bus information (needed in pommes)\n",
    "        overall_dict[col+\"_\"+year][col+\"_country\"] = \"DE\"\n",
    "        overall_dict[col+\"_\"+year][col+\"_from\"] = \"DE_bus_el\"\n",
    "        \n",
    "        \n",
    "        # Rounding: costs to 1 digit; durations to nearest integer; other params to 0 digits\n",
    "        cols_duration = [i+\"_\"+j for i in [col] for j in dur_cols]\n",
    "        cols_costs = [i+\"_\"+j for i in [col] for j in cost_cols]\n",
    "        cols_other = [i+\"_\"+j for i in [col] for j in other_cols]\n",
    "        \n",
    "        # Round durations and replace zero values\n",
    "        overall_dict[col+\"_\"+year][cols_duration] = overall_dict[col+\"_\"+year][cols_duration].round(0)\n",
    "        overall_dict[col+\"_\"+year][cols_duration] = overall_dict[col+\"_\"+year][cols_duration].replace({0:1})\n",
    "        overall_dict[col+\"_\"+year][cols_costs] = overall_dict[col+\"_\"+year][cols_costs].round(1)\n",
    "        # Replace zero cost values which might occur in rounding\n",
    "        overall_dict[col+\"_\"+year][cols_costs] = overall_dict[col+\"_\"+year][cols_costs].replace({0.0:0.01})\n",
    "        overall_dict[col+\"_\"+year][cols_other] = overall_dict[col+\"_\"+year][cols_other].round(0)\n",
    "\n",
    "        # Do some renaming\n",
    "        overall_dict[col+\"_\"+year].rename(name_dict, inplace=True)\n",
    "        overall_dict[col+\"_\"+year].drop(columns=\"sector\", inplace=True)\n",
    "        overall_dict[col+\"_\"+year].columns = overall_dict[col+\"_\"+year].columns.str.split(\n",
    "            \"_\", 1, expand=True\n",
    "        ).get_level_values(1)\n",
    "        overall_dict[col+\"_\"+year].rename(columns={np.nan: \"kind\"})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "is_executing": false,
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# Info on remaining categories is optionally stored to match it with availability time series\n",
    "# resp. to assign similar availability time series when data is missing.\n",
    "parameters_for_clustering = parameters_for_clustering.reset_index(level=2)\n",
    "\n",
    "if write_categories:\n",
    "    parameters_for_clustering[\n",
    "        parameters_for_clustering[\"year\"] == \"SQ\"\n",
    "    ].to_csv(path_folder_in+\"remaining_categories.csv\", sep=\";\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add **unit lifetime** and **efficiency** based on assumption:\n",
    "* Unit lifetime has not systematically been collected and in addition, there is only very few information from the studies evaluated.\n",
    "* Since components are mostly IT devices, a lifetime of 10 to 15 years seems realistic. The latter is used as a first proxy since even if lifetime was shorter, supposingly, there is no need for an entire reinvest, but rather a replacement of some components.\n",
    "* For efficiency, studies report only slight efficiency losses for shifting if any. Thus, it has been decided to neglect this circumstance. This is in line with not accounting for storage losses for e.g. batteries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for key in overall_dict:\n",
    "    overall_dict[key][\"unit_lifetime\"] = 15\n",
    "    overall_dict[key][\"efficiency\"] = 1.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": false
    }
   },
   "outputs": [],
   "source": [
    "# Write the parameter outputs to Excel\n",
    "if write_outputs: \n",
    "    write_multiple_sheets(overall_dict, path_folder_parameterization, filename_out+\"_overall.xlsx\") "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Write outputs\n",
    "### Write status quo data\n",
    "Write outputs needed for the power market model runs separately to csv files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Do some adjustments:\n",
    "* Simple column naming adjustment\n",
    "* Timeseries adjustments / harmonization:\n",
    "> __*NOTE*__: _2017 is used as a simulation year for the power market model._\n",
    "> * _Availability time series were originally derived for the year 2012._\n",
    "> * _Since 2017 and 2012 both started with a Sunday, no shifts of weekdays, months or seasonal patterns is necessary._\n",
    "> * _For temperature-dependent loads, average temperature data from 2017 was derived later on and insterted into the original time series to correct for the potential bias (coming from using 2012 as well as one particular temperature-measurement)._"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Do some column naming adjustment\n",
    "cases = [col + \"_SQ\" for col in quantile_cols]\n",
    "\n",
    "for case in cases:\n",
    "    cluster_ts_by_cols[case].columns = cluster_ts_by_cols[case].columns.str.split(\n",
    "        \"_\", 1, expand=True\n",
    "    ).get_level_values(1)\n",
    "    ava_pos_ts_by_cols[case].columns = ava_pos_ts_by_cols[case].columns.str.split(\n",
    "        \"_\", 1, expand=True\n",
    "    ).get_level_values(1).str[:-4]\n",
    "    ava_neg_ts_by_cols[case].columns = ava_neg_ts_by_cols[case].columns.str.split(\n",
    "        \"_\", 1, expand=True\n",
    "    ).get_level_values(1).str[:-4]\n",
    "\n",
    "    cluster_ts_by_cols[case].rename(name_dict, axis=1, inplace=True)\n",
    "    ava_pos_ts_by_cols[case].rename(name_dict, axis=1, inplace=True)\n",
    "    ava_neg_ts_by_cols[case].rename(name_dict, axis=1, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Replace time series index with 2017 one\n",
    "new_timeindex = pd.date_range(start=\"2017-01-01 00:00\", periods=8784, freq=\"H\")\n",
    "\n",
    "for case in cases:\n",
    "    cluster_ts_by_cols[case][\"new_timeindex\"] = new_timeindex\n",
    "    ava_pos_ts_by_cols[case][\"new_timeindex\"] = new_timeindex\n",
    "    ava_neg_ts_by_cols[case][\"new_timeindex\"] = new_timeindex\n",
    "    \n",
    "    cluster_ts_by_cols[case] = cluster_ts_by_cols[case].set_index(\n",
    "        \"new_timeindex\", drop=True\n",
    "    ).iloc[:8760].round(4)\n",
    "    ava_pos_ts_by_cols[case] = ava_pos_ts_by_cols[case].set_index(\n",
    "        \"new_timeindex\", drop=True\n",
    "    ).iloc[:8760].round(4)\n",
    "    ava_neg_ts_by_cols[case] = ava_neg_ts_by_cols[case].set_index(\n",
    "        \"new_timeindex\", drop=True\n",
    "    ).iloc[:8760].round(4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cluster_ts_by_cols[\"5%_SQ\"].head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write STATUS QUO data separately\n",
    "for case in cases:\n",
    "    # Node information\n",
    "    overall_dict[case].to_csv(\n",
    "        path_folder_parameterization+\"sinks_demand_response_el_\"\n",
    "        +case.split(\"_\")[0][:-1]+\".csv\"\n",
    "    )\n",
    "    # Load profile information\n",
    "    cluster_ts_by_cols[case].to_csv(\n",
    "        path_folder_parameterization+\"sinks_demand_response_el_ts_\"\n",
    "        +case.split(\"_\")[0][:-1]+\".csv\"\n",
    "    )\n",
    "    # Availability information\n",
    "    ava_pos_ts_by_cols[case].to_csv(\n",
    "        path_folder_parameterization+\"sinks_demand_response_el_ava_pos_ts_\"\n",
    "        +case.split(\"_\")[0][:-1]+\".csv\"\n",
    "    )\n",
    "    ava_neg_ts_by_cols[case].to_csv(\n",
    "        path_folder_parameterization+\"sinks_demand_response_el_ava_neg_ts_\"\n",
    "        +case.split(\"_\")[0][:-1]+\".csv\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Rearrange data to create time series per cluster\n",
    "* Extract data from `overall_dict` which contains all potential parameter information, grouped by quantile (5%, 50%, 95%) as well as year (SQ, 2020, 2025, ..., 2050).\n",
    "* Rearrange to time_series for the different parameters per demand response cluster.\n",
    "* Write the obtained results to csv."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "estimates = {}\n",
    "estimates[quantile_cols[0]] = [\n",
    "    key for key in overall_dict.keys() \n",
    "    if quantile_cols[0] in key and quantile_cols[-1] not in key\n",
    "]\n",
    "estimates[quantile_cols[1]] = [key for key in overall_dict.keys() if quantile_cols[1] in key]\n",
    "estimates[quantile_cols[-1]] = [key for key in overall_dict.keys() if quantile_cols[-1] in key]\n",
    "\n",
    "costs_columns = [\"fixed_costs\", \"specific_investments\", \"variable_costs\", \"variable_costs_shed\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for cluster in name_dict.values():\n",
    "    for scenario, estimates_list in estimates.items():\n",
    "        to_concat = []\n",
    "        for estimate in estimates_list:\n",
    "            year = estimate.split(\"_\")[1]\n",
    "            # display(overall_dict[estimate])\n",
    "            entry = overall_dict[estimate].loc[[cluster]]\n",
    "            entry[\"year\"] = year\n",
    "            to_concat.append(entry)\n",
    "        cluster_potential_parameter_data = pd.concat(to_concat)\n",
    "        cluster_potential_parameter_data.to_csv(\n",
    "            f\"{path_folder_parameterization}{cluster}_potential_parameters_real_{scenario}.csv\"\n",
    "        )\n",
    "        # Derive nominal values, whereby status quo is treated as values for 2020 (i.e. no inflation)\n",
    "        cluster_potential_parameter_data_nominal = cluster_potential_parameter_data.copy()\n",
    "        for costs_column in costs_columns:\n",
    "            cluster_potential_parameter_data_nominal.loc[\n",
    "                cluster_potential_parameter_data_nominal.year != \"SQ\", costs_column\n",
    "            ] = (\n",
    "                    cluster_potential_parameter_data.loc[\n",
    "                        cluster_potential_parameter_data_nominal.year != \"SQ\", costs_column\n",
    "                    ] \n",
    "                    * (1 + inflation_rate) ** (cluster_potential_parameter_data.loc[\n",
    "                        cluster_potential_parameter_data_nominal.year != \"SQ\", \"year\"\n",
    "                    ].astype(int) - 2020)\n",
    "                )\n",
    "        cluster_potential_parameter_data_nominal.round(4).to_csv(\n",
    "            f\"{path_folder_parameterization}{cluster}_potential_parameters_nominal_{scenario}.csv\"\n",
    "        )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bibliography\n",
    "Benz, Fabian (2019): Ermittlung von Kosten und Zeitverfügbarkeit für flexible Stromnachfragen in der Industrie in Deutschland, Freie wissenschaftliche Arbeit zur Erlangung des Grades eines Bachelor of Science am Fachgebiet Energie- und Ressourcenmanagement der TU Berlin.\n",
    "\n",
    "Gils, Hans Christian (2015): Balancing of Intermittent Renewable Power Generation by Demand Response and Thermal Energy Storage. Dissertation. Universität Stuttgart, Stuttgart.\n",
    "\n",
    "Kochems, Johannes (2020): Lastflexibilisierungspotenziale in Deutschland - Bestandsaufnahme und Entwicklungsprojektionen. Langfassung. In: IEE TU Graz (Hg.): EnInnov 2020 - 16. Symposium Energieinnovation. Energy for Future - Wege zur Klimaneutralität. Graz, 12.-14.02, https://www.tugraz.at/fileadmin/user_upload/tugrazExternal/4778f047-2e50-4e9e-b72d-e5af373f95a4/files/lf/Session_E5/553_LF_Kochems.pdf, accessed 11.05.2020.\n",
    "\n",
    "Odeh, Jonas (2019): Ermittlung von Kosten und Zeitverfügbarkeiten für die Flexibilisierung von Stromnachfragen im GHD-Sektor, Freie wissenschaftliche Arbeit zur Erlangung des Grades eines Bachelor of Science am Fachgebiet Energie- und Ressourcenmanagement der TU Berlin.\n",
    "\n",
    "Stange, Rico (2019): Ermittlung von Kosten und Zeitverfügbarkeiten einer Flexibilisierung der Stromnachfrage, Freie wissenschaftliche Arbeit zur Erlangung des Grades eines Bachelor of Science am Fachgebiet Energie- und Ressourcenmanagement der TU Berlin."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.13"
  },
  "pycharm": {
   "stem_cell": {
    "cell_type": "raw",
    "metadata": {
     "collapsed": false
    },
    "source": []
   }
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "540.797px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  },
  "toc-autonumbering": true,
  "toc-showcode": false,
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "oldHeight": 739.594,
   "position": {
    "height": "40px",
    "left": "1429.59px",
    "right": "20px",
    "top": "111.953px",
    "width": "475.719px"
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "varInspector_section_display": "none",
   "window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}