scenario_setup.py

# Script to run calculations to determine grid expansion costs

import os
import random
import logging
import pickle

from copy import deepcopy
import numpy as np
import pandas as pd

from data_preparation.prepare_grids import \
    distribute_energy_resources
from edisgo.opf.timeseries_reduction import get_steps_reinforcement
from edisgo.edisgo import import_edisgo_from_pickle
from Script_prepare_grids_for_optimization import get_downstream_nodes_matrix_iterative

from edisgo import EDisGo
from edisgo.network.results import Results
from edisgo.network.timeseries import TimeSeries
from edisgo.tools.complexity_reduction import remove_1m_lines_from_edisgo

from data_preparation import generator_scenario
from results_visualisation import load_costs

USER_BASEPATH = os.path.dirname(__file__)

logger = logging.getLogger(__name__)
logger.setLevel(logging.DEBUG)


def scenario_input_data():
    """
    Returns dictionary with relevant scenario data.

    Returns
    --------
    dict
        Dictionary with relevant scenario data. Should be changed manually to adapt
        scenario.
    """
    return {
        "PV_installed_capacities": 0.008, # MW
        "PV_gamma": (190.88895608306424, -46.093268849468814, 0.2853110801380861),
        "BESS_relative_capacity": 1.0, # kWh/kW_PV
        "BESS_relative_power": 0.6, # kW/kWh
        "BESS_charging_efficiency": 1.0,
        "BESS_discharging_efficiency": 1.0,
        "HP_installed_capacities": 0.002, # MW
        "HP_share_air_source": 0.8,
        "Annual_heat_demand_house": 140*45*1e-3, # area in m^2 times heat demand in
                                                 # kWh/m^2 p.a. converted into MWh p.a
        "HP_size_TES": 0.016, # MWh
        "HP_gamma": (5.433654509183175, -0.770275842540082, 2.543186523886985),
    }


def get_seed_with_min_rmsd(results, scenario, grid_id):
    results_tmp = results.loc[(results.scenario == scenario) & (results.grid_id==grid_id)]
    results_mean = results_tmp.groupby("penetration").mean()
    rmsd = pd.DataFrame(index=results_tmp.seed.unique(), columns=["rmsd"])
    for seed in results_tmp.seed.unique():
        rmsd.loc[seed, "rmsd"] = np.sqrt(
            (results_tmp.loc[results_tmp.seed==seed, "costs"]-results_mean.loc[results_tmp.loc[results_tmp.seed==seed, "penetration"], "costs"].values).
            apply(np.square).sum())
    return rmsd[rmsd==rmsd.min()].dropna(how="all").index[0]


def save_added_components(results_dir, resources, edisgo_obj, edisgo_orig):
    if "RES" in resources:
        added_res = edisgo_obj.topology.generators_df.loc[~edisgo_obj.topology.generators_df.index.isin(
            edisgo_orig.topology.generators_df.index)]
        added_res.to_csv(os.path.join(results_dir, "added_generators.csv"))
    if "BESS" in resources:
        added_bess = edisgo_obj.topology.storage_units_df.loc[~edisgo_obj.topology.storage_units_df.index.isin(
            edisgo_orig.topology.storage_units_df.index)]
        added_bess.to_csv(os.path.join(results_dir, "added_storage_units.csv"))
    if "HP" in resources:
        added_hps = edisgo_obj.topology.loads_df.loc[~edisgo_obj.topology.loads_df.index.isin(
            edisgo_orig.topology.loads_df.index) & (edisgo_obj.topology.loads_df.type == "heat_pump")]
        added_hps.to_csv(os.path.join(results_dir, "added_heat_pumps.csv"))
    if "EV" in resources:
        added_evs = edisgo_obj.topology.loads_df.loc[~edisgo_obj.topology.loads_df.index.isin(
            edisgo_orig.topology.loads_df.index) & (edisgo_obj.topology.loads_df.type == "charging_point")]
        added_evs.to_csv(os.path.join(results_dir, "added_charging_points.csv"))


def setup_timeindex(grid_id, resources, run_full, ts_mode=None, max_day=None):
    # use timesteps of reference run
    timesteps = pd.date_range("2011-1-1", periods=0, freq="H")
    simulated_days = []
    for resource in resources:
        if ts_mode == "debug":
            if max_day is not None:
                logger.warning(
                    "Max_day parameter defined in debug mode. In debug mode, only one "
                    "day per resource will be simulated. Please set debug to False if "
                    "more days should be simulated.")
            if resource == "RES":
                time_extreme = pd.to_datetime("2011-03-07 13:00:00")
            elif resource == "HP":
                time_extreme = pd.to_datetime("2011-02-23 05:00:00")
            elif resource == "EV":
                time_extreme = pd.to_datetime("2011-08-23 18:00:00")
            elif resource == "BESS":
                continue
            else:
                raise ValueError(f"The resource you entered is not valid: {resource}.")
            timesteps = timesteps.union(pd.date_range(
                time_extreme.date(), periods=24, freq="H"))
        else:
            if resource == "BESS":
                resource = "RES_BESS"
            scenario = f"{resource}_reference"
            try:
                relevant_timesteps_tmp = \
                    pd.read_csv(os.path.join(run_full["result_dir"], run_full["run"], scenario,
                                             str(grid_id), str(run_full["penetration"]),
                                             str(run_full["seed"]),
                                             "relevant_timesteps.csv"),
                                index_col=0, parse_dates=True)
            except FileNotFoundError:
                raise FileNotFoundError("You have to provide a directory in which the "
                                        "relevant time steps are saved. Run the "
                                        "simulation in full_ts mode to determine and "
                                        "save the relevant steps.")
            if max_day is None:
                max_day = len(relevant_timesteps_tmp.index)
            for timestep in relevant_timesteps_tmp.index:
                if not timestep.date() in simulated_days:
                    timesteps = timesteps.union(pd.date_range(timestep.date(), periods=24, freq="H"))
                    simulated_days.append(timestep.date())
                # if maximum number of days is reached, break
                if len(simulated_days) >= max_day:
                    break
    return timesteps


def setup_original_grid(grids_path, grid_id, scenario_dict):
    grid_path = os.path.join(grids_path, str(grid_id))
    # try to load original grid
    try:
        edisgo_obj_orig = \
            import_edisgo_from_pickle("edisgo_object.pkl", grid_path)
    # if grid does not exist yet, create and save original grid
    except FileNotFoundError:
        os.makedirs(grid_path, exist_ok=True)
        edisgo_obj_orig = setup_topology_status_quo_ding0(
            scenario_dict["grids_orig_path"], grid_id)
        edisgo_obj_orig.save_edisgo_to_pickle(grid_path, "edisgo_object.pkl")
        downstream_matrix = \
            get_downstream_nodes_matrix_iterative(edisgo_obj_orig.topology)
        pickle.dump(downstream_matrix,
                    open(os.path.join(grid_path, "downstream_matrix.pkl"), "wb"))
    if not scenario_dict["ts_mode"] == "full":
        edisgo_obj_orig.timeseries.timeindex = setup_timeindex(
            grid_id=grid_id,
            resources=scenario_dict["resources"],
            run_full=scenario_dict["run_full"],
            max_day=scenario_dict["max_day"],
            ts_mode=scenario_dict["ts_mode"],
        )
        # adapt timeseries
        attributes = TimeSeries()._attributes
        for attr in attributes:
            if not getattr(edisgo_obj_orig.timeseries, attr).empty:
                setattr(
                    edisgo_obj_orig.timeseries, attr,
                    getattr(edisgo_obj_orig.timeseries, attr).loc[edisgo_obj_orig.timeseries.timeindex]
                )
    return edisgo_obj_orig


def setup_topology_status_quo_ding0(grids_path, grid_id):
    """
    Reinforces the grids to be stable in the status quo.

    Generator capacities (wind, solar, biomass) are updated to better reflect installed
    capacities in 2018 (generator capacities in ding0 grids are from 2015). Afterwards,
    grids are reinforced.

    Parameters
    -----------
    grids_path : str
        Path to directory containing original ding0 grids.
    grid_id : int or str
        ID of ding0 grid to set up.

    Returns
    ---------
    EDisGo object
        EDisGo object with status quo grid topology and worst case time series.

    """

    dingo_grid_path = os.path.join(grids_path, str(grid_id))
    edisgo_obj = EDisGo(
        ding0_grid=dingo_grid_path,
        config_path=os.path.join(USER_BASEPATH, "data_preparation", "data", "config_data")
    )

    # install 2018 generator capacities
    logger.debug("Installed capacities 2015: {}".format(
        edisgo_obj.topology.generators_df.loc[
            :, ["p_nom", "type"]].groupby(["type"]).sum().loc[
            :, "p_nom"]
    ))
    target_capacities = generator_scenario.set_up_status_quo_2018(
        edisgo_obj
    )
    edisgo_obj.import_generators(
        generator_scenario="nep2035",
        p_target=target_capacities,
        remove_decommissioned=False,
        update_existing=False
    )
    logger.debug("Installed capacities 2018: {}".format(
        edisgo_obj.topology.generators_df.loc[
            :, ["p_nom", "type"]].groupby(["type"]).sum().loc[
            :, "p_nom"]
    ))

    # remove 1 meter lines
    edisgo_obj = remove_1m_lines_from_edisgo(edisgo_obj)

    # set up worst case time series Todo: change to timeseries data
    edisgo_obj.set_time_series_worst_case_analysis()

    # Reinforce ding0 grid to obtain a stable status quo grid
    logger.debug("Conduct grid reinforcement to obtain stable status quo grid.")
    edisgo_obj.reinforce()

    # Clear results
    edisgo_obj.results = Results(edisgo_obj)

    # Add timeseries
    edisgo_obj.timeseries = TimeSeries()
    edisgo_obj.timeseries.timeindex = \
        pd.date_range("1/1/2011 00:00", periods=8760, freq="H")
    edisgo_obj.set_time_series_active_power_predefined(
        fluctuating_generators_ts="oedb",
        dispatchable_generators_ts="full_capacity",
        conventional_loads_ts="demandlib",
    )
    edisgo_obj.set_time_series_reactive_power_control()

    # Check integrity of edisgo object
    edisgo_obj.check_integrity()

    # Reinforce with timeseries
    edisgo_obj, _ = run_reduced_reinforcement(edisgo_obj)
    # Clear results
    edisgo_obj.results = Results(edisgo_obj)

    logger.debug("Status quo grids ready.")

    return edisgo_obj


def run_calculation(grids_path, results_path, scenario_dict, seed, penetration,
                    grid_id, edisgo_obj_orig=None):
    # if run was already solved, skip
    grid_results_path = os.path.join(results_path, scenario_dict["run"],
                                     scenario_dict["scenario"], str(grid_id))
    os.makedirs(os.path.join(grid_results_path, str(penetration), str(seed)),
                exist_ok=True)
    if len(os.listdir(os.path.join(grid_results_path, str(penetration), str(seed))))>=2:
        return pd.Series()
    if edisgo_obj_orig is None:
        edisgo_obj_orig = setup_original_grid(grids_path, grid_id, scenario_dict)
    # Set new seed
    random.seed(seed)
    # Use seed and distribute DER
    edisgo_obj = distribute_energy_resources(
        edisgo_obj=deepcopy(edisgo_obj_orig),
        scenario_dict=scenario_dict,
        seed=seed,
        penetration=penetration,
        resources=scenario_dict["resources"],
        strategy=scenario_dict["strategy"],
        ev_dir=os.path.join(scenario_dict["ev_base_dir"], str(grid_id)),
        sizing_mode=scenario_dict["sizing_mode"],
        matrix_path=grids_path
    )
    edisgo_obj.check_integrity()
    # edisgo_obj.save_edisgo_to_pickle(grid_path, "edisgo_object_EV_ref.pkl")
    # Analyse grid and save results
    edisgo_obj, relevant_timesteps = run_reduced_reinforcement(
        edisgo_obj=edisgo_obj,
        num_steps_current=scenario_dict["num_ts_loading"],
        num_steps_voltage=scenario_dict["num_ts_voltage"],
    )
    # save results
    pd.DataFrame(index=relevant_timesteps).to_csv(
        os.path.join(grid_results_path, str(penetration), str(seed),
                     "relevant_timesteps.csv"))
    save_added_components(
        results_dir=os.path.join(grid_results_path, str(penetration), str(seed)),
        resources=scenario_dict["resources"],
        edisgo_obj=edisgo_obj,
        edisgo_orig=edisgo_obj_orig,
    )
    # if "BESS" in scenario_dict["resources"]:
    #     os.makedirs(os.path.join(grid_results_path, str(penetration), str(seed)),
    #                 exist_ok=True)
    #     # edisgo_obj.timeseries.storage_units_active_power.to_csv(
    #     #     os.path.join(grid_results_path, str(penetration), str(seed),
    #     #                  "storage_units_active_power.csv"))
    # extract grid expansion costs
    if not edisgo_obj.results.grid_expansion_costs.empty:
        grid_expansion_costs = \
            edisgo_obj.results.grid_expansion_costs.total_costs.sum()
        edisgo_obj.results.grid_expansion_costs.to_csv(
            os.path.join(grid_results_path, str(penetration), str(seed),
                         "grid_expansion_costs.csv")
        )
    else:
        grid_expansion_costs = 0.0
    costs = pd.Series({
                        "grid_id": grid_id,
                        "seed": seed,
                        "penetration": penetration,
                        "costs": grid_expansion_costs
                    })
    return costs


def run_reduced_reinforcement(edisgo_obj, num_steps_current=None,
                              num_steps_voltage=None):
    # Todo: implement alternative when no timesteps converge
    timesteps_not_converged = [0]
    reinforcement_timesteps = pd.DatetimeIndex([])
    while len(timesteps_not_converged) > 0:
        timesteps_not_converged = edisgo_obj.analyze(raise_not_converged=False)
        if len(timesteps_not_converged) == len(edisgo_obj.timeseries.timeindex):
            print("WARNING: None of the timesteps converged in power flow.")
            return edisgo_obj, timesteps_not_converged
        # Reinforce grid and save results
        relevant_steps = get_steps_reinforcement(
            edisgo_obj=edisgo_obj,
            num_steps_current=num_steps_current,
            num_steps_voltage=num_steps_voltage
        )
        # Todo: combined_analysis=True? Which factors for reinforcement?
        if len(relevant_steps) > 0:
            edisgo_obj.reinforce(timesteps_pfa=relevant_steps)
            reinforcement_timesteps = reinforcement_timesteps.append(relevant_steps)
        else:
            logger.info("No reinforcement necessary.")
        # Todo: save timeseries or not? Maybe only save timeseries of der
    return edisgo_obj, list(set(reinforcement_timesteps))


if __name__ == "__main__":
    res_dir = r"H:\cost_functions\results"
    runs = ["Powertech_2022"]
    scenarios = ["HP_reference", "EV_reference", "RES_reference", "RES_BESS_reference", "RES_BESS_EV_HP_reference"]
    grid_ids = [176, 177, 1056, 1690, 1811, 2534]
    penetrations = [
        0,
        0.1,
        0.2,
        0.3,
        0.4,
        0.5,
        0.6,
        0.7,
        0.8,
        0.9,
        1.0
    ]
    seeds = [
        544,
        295,
        453,
        558,
        317,
        599,
        62,
        530,
        802,
        708]
    results = load_costs(res_dir, runs, scenarios, grid_ids, penetrations, seeds)
    representative_seeds = pd.DataFrame(index=grid_ids, columns=scenarios)
    for scenario in scenarios:
        for grid_id in grid_ids:
            representative_seeds.loc[grid_id, scenario] = \
                get_seed_with_min_rmsd(results, scenario, grid_id)
    representative_seeds.to_csv(os.path.join(res_dir, runs[0], "representative_seeds.csv"))
    print("Success")