diff --git a/edisgo/flex_opt/check_tech_constraints.py b/edisgo/flex_opt/check_tech_constraints.py
index 0bbc91378..ccba6def0 100644
--- a/edisgo/flex_opt/check_tech_constraints.py
+++ b/edisgo/flex_opt/check_tech_constraints.py
@@ -1,8 +1,9 @@
-import pandas as pd
+import itertools
 import logging
 from math import sqrt
+
 import numpy as np
-import itertools
+import pandas as pd
 
 from edisgo.network.grids import LVGrid, MVGrid
 
@@ -136,25 +137,29 @@ def lines_allowed_load(edisgo_obj, voltage_level):
         )
 
     i_lines_allowed_per_case = {}
-    i_lines_allowed_per_case["feedin_case"] = (
+    i_lines_allowed_per_case["feed-in_case"] = (
         lines_df.s_nom
         / sqrt(3)
         / nominal_voltage
         * edisgo_obj.config["grid_expansion_load_factors"][
-            "{}_feedin_case_line".format(voltage_level)
+            "{}_feed-in_case_line".format(voltage_level)
         ]
     )
 
     # adapt i_lines_allowed for radial feeders
     buses_in_cycles = list(
-        set(itertools.chain.from_iterable(edisgo_obj.topology.rings)))
+        set(itertools.chain.from_iterable(edisgo_obj.topology.rings))
+    )
 
     # Find lines in cycles
     lines_in_cycles = list(
-        lines_df.loc[lines_df[[
-            'bus0', 'bus1']].isin(buses_in_cycles).all(axis=1)].index.values)
+        lines_df.loc[
+            lines_df[["bus0", "bus1"]].isin(buses_in_cycles).all(axis=1)
+        ].index.values
+    )
     lines_radial_feeders = list(
-        lines_df.loc[~lines_df.index.isin(lines_in_cycles)].index.values)
+        lines_df.loc[~lines_df.index.isin(lines_in_cycles)].index.values
+    )
 
     # lines in cycles have to be n-1 secure
     i_lines_allowed_per_case["load_case"] = (
@@ -167,9 +172,9 @@ def lines_allowed_load(edisgo_obj, voltage_level):
     )
 
     # lines in radial feeders are not n-1 secure anyways
-    i_lines_allowed_per_case["load_case"] = \
-        i_lines_allowed_per_case["load_case"].append(
-            lines_df.loc[lines_radial_feeders].s_nom / sqrt(3) / nominal_voltage)
+    i_lines_allowed_per_case["load_case"] = i_lines_allowed_per_case[
+        "load_case"
+    ].append(lines_df.loc[lines_radial_feeders].s_nom / sqrt(3) / nominal_voltage)
 
     i_lines_allowed = edisgo_obj.timeseries.timesteps_load_feedin_case.loc[
         edisgo_obj.results.i_res.index
@@ -198,8 +203,9 @@ def lines_relative_load(edisgo_obj, lines_allowed_load):
 
     """
     # get line load from power flow analysis
-    i_lines_pfa = edisgo_obj.results.i_res.loc[lines_allowed_load.index,
-                                               lines_allowed_load.columns]
+    i_lines_pfa = edisgo_obj.results.i_res.loc[
+        lines_allowed_load.index, lines_allowed_load.columns
+    ]
 
     return i_lines_pfa / lines_allowed_load
 
@@ -241,9 +247,7 @@ def _line_load(edisgo_obj, voltage_level):
     # calculate relative line load and keep maximum over-load of each line
     relative_i_res = lines_relative_load(edisgo_obj, i_lines_allowed)
 
-    crit_lines_relative_load = (
-        relative_i_res[relative_i_res > 1].max().dropna()
-    )
+    crit_lines_relative_load = relative_i_res[relative_i_res > 1].max().dropna()
 
     if len(crit_lines_relative_load) > 0:
         crit_lines = pd.concat(
@@ -253,7 +257,7 @@ def _line_load(edisgo_obj, voltage_level):
             ],
             axis=1,
             keys=["max_rel_overload", "time_index"],
-            sort=True
+            sort=True,
         )
         crit_lines.loc[:, "voltage_level"] = voltage_level
     else:
@@ -326,9 +330,7 @@ def mv_lv_station_load(edisgo_obj):
 
     crit_stations = pd.DataFrame()
     for lv_grid in edisgo_obj.topology.mv_grid.lv_grids:
-        crit_stations = crit_stations.append(
-            _station_load(edisgo_obj, lv_grid)
-        )
+        crit_stations = crit_stations.append(_station_load(edisgo_obj, lv_grid))
     if not crit_stations.empty:
         logger.debug(
             "==> {} MV/LV station(s) has/have load issues.".format(
@@ -366,9 +368,9 @@ def _station_load(edisgo_obj, grid):
     if isinstance(grid, LVGrid):
         voltage_level = "lv"
         transformers_df = grid.transformers_df
-        s_station_pfa = edisgo_obj.results.s_res.loc[
-            :, transformers_df.index
-        ].sum(axis=1)
+        s_station_pfa = edisgo_obj.results.s_res.loc[:, transformers_df.index].sum(
+            axis=1
+        )
     elif isinstance(grid, MVGrid):
         voltage_level = "mv"
         transformers_df = edisgo_obj.topology.transformers_hvmv_df
@@ -377,7 +379,8 @@ def _station_load(edisgo_obj, grid):
         if not any(mv_lines.isin(edisgo_obj.results.i_res.columns)):
             raise ValueError(
                 "MV was not included in power flow analysis, wherefore load "
-                "of HV/MV station cannot be calculated.")
+                "of HV/MV station cannot be calculated."
+            )
         s_station_pfa = np.hypot(
             edisgo_obj.results.pfa_slack.p,
             edisgo_obj.results.pfa_slack.q,
@@ -534,9 +537,7 @@ def lv_voltage_deviation(edisgo_obj, mode=None, voltage_levels="mv_lv"):
     crit_buses = {}
 
     if voltage_levels == "mv_lv":
-        v_limits_upper, v_limits_lower = _mv_allowed_voltage_limits(
-            edisgo_obj, "mv_lv"
-        )
+        v_limits_upper, v_limits_lower = _mv_allowed_voltage_limits(edisgo_obj, "mv_lv")
     elif not "lv" == voltage_levels:
         raise ValueError(
             "{} is not a valid option for input variable 'voltage_levels' in "
@@ -573,15 +574,11 @@ def lv_voltage_deviation(edisgo_obj, mode=None, voltage_levels="mv_lv"):
     if crit_buses:
         if mode == "stations":
             logger.debug(
-                "==> {} LV station(s) has/have voltage issues.".format(
-                    len(crit_buses)
-                )
+                "==> {} LV station(s) has/have voltage issues.".format(len(crit_buses))
             )
         else:
             logger.debug(
-                "==> {} LV topology(s) has/have voltage issues.".format(
-                    len(crit_buses)
-                )
+                "==> {} LV topology(s) has/have voltage issues.".format(len(crit_buses))
             )
     else:
         if mode == "stations":
@@ -627,9 +624,9 @@ def _mv_allowed_voltage_limits(edisgo_obj, voltage_levels):
 
     # get config values for lower voltage limit in feed-in case and upper
     # voltage limit in load case
-    v_allowed_per_case["feedin_case_lower"] = edisgo_obj.config[
+    v_allowed_per_case["feed-in_case_lower"] = edisgo_obj.config[
         "grid_expansion_allowed_voltage_deviations"
-    ]["feedin_case_lower"]
+    ]["feed-in_case_lower"]
     v_allowed_per_case["load_case_upper"] = edisgo_obj.config[
         "grid_expansion_allowed_voltage_deviations"
     ]["load_case_upper"]
@@ -639,17 +636,17 @@ def _mv_allowed_voltage_limits(edisgo_obj, voltage_levels):
     offset = edisgo_obj.config["grid_expansion_allowed_voltage_deviations"][
         "hv_mv_trafo_offset"
     ]
-    control_deviation = edisgo_obj.config[
-        "grid_expansion_allowed_voltage_deviations"
-    ]["hv_mv_trafo_control_deviation"]
+    control_deviation = edisgo_obj.config["grid_expansion_allowed_voltage_deviations"][
+        "hv_mv_trafo_control_deviation"
+    ]
 
     if voltage_levels == "mv_lv" or voltage_levels == "mv":
-        v_allowed_per_case["feedin_case_upper"] = (
+        v_allowed_per_case["feed-in_case_upper"] = (
             1
             + offset
             + control_deviation
             + edisgo_obj.config["grid_expansion_allowed_voltage_deviations"][
-                "{}_feedin_case_max_v_deviation".format(voltage_levels)
+                "{}_feed-in_case_max_v_deviation".format(voltage_levels)
             ]
         )
         v_allowed_per_case["load_case_lower"] = (
@@ -715,17 +712,15 @@ def _lv_allowed_voltage_limits(edisgo_obj, lv_grid, mode):
         config_string = "mv_lv_station"
     else:
         # reference voltage is voltage at stations' secondary side
-        voltage_base = edisgo_obj.results.v_res.loc[
-            :, lv_grid.station.index.values[0]
-        ]
+        voltage_base = edisgo_obj.results.v_res.loc[:, lv_grid.station.index.values[0]]
         config_string = "lv"
 
     # calculate upper voltage limit in feed-in case and lower voltage limit in
     # load case
-    v_allowed_per_case["feedin_case_upper"] = (
+    v_allowed_per_case["feed-in_case_upper"] = (
         voltage_base
         + edisgo_obj.config["grid_expansion_allowed_voltage_deviations"][
-            "{}_feedin_case_max_v_deviation".format(config_string)
+            "{}_feed-in_case_max_v_deviation".format(config_string)
         ]
     )
     v_allowed_per_case["load_case_lower"] = (
@@ -736,9 +731,9 @@ def _lv_allowed_voltage_limits(edisgo_obj, lv_grid, mode):
     )
 
     timeindex = voltage_base.index
-    v_allowed_per_case["feedin_case_lower"] = pd.Series(
+    v_allowed_per_case["feed-in_case_lower"] = pd.Series(
         edisgo_obj.config["grid_expansion_allowed_voltage_deviations"][
-            "feedin_case_lower"
+            "feed-in_case_lower"
         ],
         index=timeindex,
     )
@@ -755,12 +750,8 @@ def _lv_allowed_voltage_limits(edisgo_obj, lv_grid, mode):
     load_feedin_case = edisgo_obj.timeseries.timesteps_load_feedin_case
     for t in timeindex:
         case = load_feedin_case.loc[t]
-        v_limits_upper.append(
-            v_allowed_per_case["{}_upper".format(case)].loc[t]
-        )
-        v_limits_lower.append(
-            v_allowed_per_case["{}_lower".format(case)].loc[t]
-        )
+        v_limits_upper.append(v_allowed_per_case["{}_upper".format(case)].loc[t])
+        v_limits_lower.append(v_allowed_per_case["{}_lower".format(case)].loc[t])
     v_limits_upper = pd.Series(v_limits_upper, index=timeindex)
     v_limits_lower = pd.Series(v_limits_lower, index=timeindex)
 
@@ -820,23 +811,19 @@ def voltage_diff(edisgo_obj, buses, v_dev_allowed_upper, v_dev_allowed_lower):
         (v_dev_allowed_lower.loc[v_mag_pu_pfa.index]).values,
         (v_mag_pu_pfa.shape[1], 1),
     )
-    overvoltage = v_mag_pu_pfa.T[
-        v_mag_pu_pfa.T > v_dev_allowed_upper_format
-    ].dropna(how="all")
-    undervoltage = v_mag_pu_pfa.T[
-        v_mag_pu_pfa.T < v_dev_allowed_lower_format
-    ].dropna(how="all")
+    overvoltage = v_mag_pu_pfa.T[v_mag_pu_pfa.T > v_dev_allowed_upper_format].dropna(
+        how="all"
+    )
+    undervoltage = v_mag_pu_pfa.T[v_mag_pu_pfa.T < v_dev_allowed_lower_format].dropna(
+        how="all"
+    )
     # sort buses with under- and overvoltage issues in a way that
     # worst case is saved
     buses_both = v_mag_pu_pfa[
         overvoltage[overvoltage.index.isin(undervoltage.index)].index
     ]
-    voltage_diff_ov = (
-        buses_both.T - v_dev_allowed_upper.loc[v_mag_pu_pfa.index].values
-    )
-    voltage_diff_uv = (
-        -buses_both.T + v_dev_allowed_lower.loc[v_mag_pu_pfa.index].values
-    )
+    voltage_diff_ov = buses_both.T - v_dev_allowed_upper.loc[v_mag_pu_pfa.index].values
+    voltage_diff_uv = -buses_both.T + v_dev_allowed_lower.loc[v_mag_pu_pfa.index].values
     voltage_diff_ov = voltage_diff_ov.loc[
         voltage_diff_ov.max(axis=1) > voltage_diff_uv.max(axis=1)
     ]
@@ -914,18 +901,12 @@ def _append_crit_buses(df):
 
     # append to crit buses dataframe
     if not voltage_diff_ov.empty:
-        crit_buses_grid = crit_buses_grid.append(
-            _append_crit_buses(voltage_diff_ov)
-        )
+        crit_buses_grid = crit_buses_grid.append(_append_crit_buses(voltage_diff_ov))
     if not voltage_diff_uv.empty:
-        crit_buses_grid = crit_buses_grid.append(
-            _append_crit_buses(voltage_diff_uv)
-        )
+        crit_buses_grid = crit_buses_grid.append(_append_crit_buses(voltage_diff_uv))
 
     if not crit_buses_grid.empty:
-        crit_buses_grid.sort_values(
-            by=["v_diff_max"], ascending=False, inplace=True
-        )
+        crit_buses_grid.sort_values(by=["v_diff_max"], ascending=False, inplace=True)
 
     return crit_buses_grid
 
diff --git a/edisgo/network/electromobility.py b/edisgo/network/electromobility.py
index bf037652d..7ec6a7ca1 100644
--- a/edisgo/network/electromobility.py
+++ b/edisgo/network/electromobility.py
@@ -1,15 +1,14 @@
-import os
 import logging
-import pandas as pd
-import numpy as np
-import geopandas as gpd
-
+import os
 from zipfile import ZipFile
+
+import geopandas as gpd
+import numpy as np
+import pandas as pd
 from sklearn import preprocessing
 
 from edisgo.network.components import PotentialChargingParks
 
-
 logger = logging.getLogger("edisgo")
 
 COLUMNS = {
@@ -23,12 +22,7 @@
         "park_start",
         "park_end",
     ],
-    "grid_connections_gdf": [
-        "id",
-        "use_case",
-        "user_centric_weight",
-        "geometry"
-    ],
+    "grid_connections_gdf": ["id", "use_case", "user_centric_weight", "geometry"],
     "simbev_config_df": ["value"],
     "potential_charging_parks_df": [
         "lv_grid_id",
@@ -145,8 +139,7 @@ def integrated_charging_parks_df(self):
         try:
             return self._integrated_charging_parks_df
         except Exception:
-            return pd.DataFrame(
-                columns=COLUMNS["integrated_charging_parks_df"])
+            return pd.DataFrame(columns=COLUMNS["integrated_charging_parks_df"])
 
     @integrated_charging_parks_df.setter
     def integrated_charging_parks_df(self, df):
@@ -236,8 +229,7 @@ def to_csv(self, directory):
                 path = os.path.join(directory, file)
                 df.to_csv(path)
 
-    def from_csv(
-            self, data_path, edisgo_obj, from_zip_archive=False):
+    def from_csv(self, data_path, edisgo_obj, from_zip_archive=False):
         """
         Restores electromobility from csv files.
 
@@ -274,33 +266,32 @@ def from_csv(
             if from_zip_archive:
                 # open zip file to make it readable for pandas
                 with zip.open(file) as f:
-                    df = pd.read_csv(
-                        f, index_col=0)
+                    df = pd.read_csv(f, index_col=0)
             else:
                 path = os.path.join(data_path, file)
-                df = pd.read_csv(
-                    path, index_col=0)
+                df = pd.read_csv(path, index_col=0)
 
             if attr == "grid_connections_gdf":
                 epsg = edisgo_obj.topology.grid_district["srid"]
 
-                df = df.assign(
-                    geometry=gpd.GeoSeries.from_wkt(df["geometry"]))
+                df = df.assign(geometry=gpd.GeoSeries.from_wkt(df["geometry"]))
 
                 try:
                     df = gpd.GeoDataFrame(
-                        df, geometry="geometry", crs={"init": f"epsg:{epsg}"})
+                        df, geometry="geometry", crs={"init": f"epsg:{epsg}"}
+                    )
 
                 except Exception as _:
                     logging.warning(
                         f"Grid connections could not be loaded with "
-                        f"EPSG {epsg}. Trying with EPSG 4326 as fallback.")
+                        f"EPSG {epsg}. Trying with EPSG 4326 as fallback."
+                    )
 
                     df = gpd.GeoDataFrame(
-                        df, geometry="geometry", crs={"init": "epsg:4326"})
+                        df, geometry="geometry", crs={"init": "epsg:4326"}
+                    )
 
-            setattr(
-                self, attr, df)
+            setattr(self, attr, df)
 
         if from_zip_archive:
             # make sure to destroy ZipFile Class to close any open connections
@@ -328,34 +319,31 @@ def _potential_charging_parks_df(self):
             potential_charging_parks = list(self.potential_charging_parks)
 
             potential_charging_parks_df.lv_grid_id = [
-                _.nearest_substation["lv_grid_id"]
-                for _ in potential_charging_parks
+                _.nearest_substation["lv_grid_id"] for _ in potential_charging_parks
             ]
 
             potential_charging_parks_df.distance_to_nearest_substation = [
-                _.nearest_substation["distance"]
-                for _ in potential_charging_parks
+                _.nearest_substation["distance"] for _ in potential_charging_parks
             ]
 
             min_max_scaler = preprocessing.MinMaxScaler()
 
-            potential_charging_parks_df.distance_weight = (
-                1 - min_max_scaler.fit_transform(
-                    potential_charging_parks_df.
-                    distance_to_nearest_substation.values.reshape(-1, 1)
+            potential_charging_parks_df.distance_weight = 1 - min_max_scaler.fit_transform(
+                potential_charging_parks_df.distance_to_nearest_substation.values.reshape(
+                    -1, 1
                 )
             )
 
             potential_charging_parks_df.charging_point_capacity = [
-                _.designated_charging_point_capacity
-                for _ in potential_charging_parks
+                _.designated_charging_point_capacity for _ in potential_charging_parks
             ]
 
             potential_charging_parks_df.charging_point_weight = (
                 1
                 - min_max_scaler.fit_transform(
-                    potential_charging_parks_df.charging_point_capacity.
-                    values.reshape(-1, 1)
+                    potential_charging_parks_df.charging_point_capacity.values.reshape(
+                        -1, 1
+                    )
                 )
             )
 
@@ -389,16 +377,17 @@ def _get_matching_dict_of_attributes_and_file_names():
     return emob_dict
 
 
-def get_energy_bands_for_optimization(edisgo_obj, use_case, t_max=372*4*24):
+def get_energy_bands_for_optimization(edisgo_obj, use_case, t_max=372 * 4 * 24):
     """
     Method to extract flexibility bands for linear optimisation.
     """
+
     def _shorten_and_set_index(band):
         """
         Method to adjust bands to timeindex of edisgo_obj
         #Todo: change such that first day is replaced by (365+1)th day
         """
-        band = band.iloc[:len(edisgo_obj.timeseries.timeindex)]
+        band = band.iloc[: len(edisgo_obj.timeseries.timeindex)]
         band.index = edisgo_obj.timeseries.timeindex
         return band
 
@@ -408,12 +397,9 @@ def _shorten_and_set_index(band):
     ]
     # set up bands
     tmp_idx = range(t_max)
-    upper_power = pd.DataFrame(index=tmp_idx,
-                               columns=cps.index, data=0)
-    upper_energy = pd.DataFrame(index=tmp_idx,
-                                columns=cps.index, data=0)
-    lower_energy = pd.DataFrame(index=tmp_idx,
-                                columns=cps.index, data=0)
+    upper_power = pd.DataFrame(index=tmp_idx, columns=cps.index, data=0)
+    upper_energy = pd.DataFrame(index=tmp_idx, columns=cps.index, data=0)
+    lower_energy = pd.DataFrame(index=tmp_idx, columns=cps.index, data=0)
     hourly_steps = 60 / edisgo_obj.electromobility.stepsize
     for cp in cps.index:
         # get index of charging park used in charging processes
@@ -422,7 +408,9 @@ def _shorten_and_set_index(band):
         ].index
         # get relevant charging processes
         charging_processes = edisgo_obj.electromobility.charging_processes_df.loc[
-            edisgo_obj.electromobility.charging_processes_df.charging_park_id.isin(charging_park_id)
+            edisgo_obj.electromobility.charging_processes_df.charging_park_id.isin(
+                charging_park_id
+            )
         ]
         # iterate through charging processes and fill matrices
         for idx, charging_process in charging_processes.iterrows():
@@ -430,37 +418,85 @@ def _shorten_and_set_index(band):
             if charging_process.park_end == t_max:
                 continue
             # get brutto charging capacity
-            brutto_charging_capacity = charging_process.netto_charging_capacity/\
-                                       edisgo_obj.electromobility.eta_charging_points
+            brutto_charging_capacity = (
+                charging_process.netto_charging_capacity
+                / edisgo_obj.electromobility.eta_charging_points
+            )
             # charging power
-            upper_power.loc[charging_process.park_start-1:charging_process.park_end-1, cp] += \
-                brutto_charging_capacity
+            upper_power.loc[
+                charging_process.park_start - 1 : charging_process.park_end - 1, cp
+            ] += brutto_charging_capacity
             # energy bands
-            charging_time = \
-                charging_process.chargingdemand / charging_process.netto_charging_capacity * \
-                hourly_steps
-            if charging_time - (charging_process.park_end - charging_process.park_start + 1) > 1e-6:
-                raise ValueError("Charging demand cannot be fulfilled for charging process {}. "
-                                 "Please check.".format(idx))
+            charging_time = (
+                charging_process.chargingdemand
+                / charging_process.netto_charging_capacity
+                * hourly_steps
+            )
+            if (
+                charging_time
+                - (charging_process.park_end - charging_process.park_start + 1)
+                > 1e-6
+            ):
+                raise ValueError(
+                    "Charging demand cannot be fulfilled for charging process {}. "
+                    "Please check.".format(idx)
+                )
             full_charging_steps = int(charging_time)
-            part_time_step = charging_time-full_charging_steps
+            part_time_step = charging_time - full_charging_steps
             # lower band
-            lower_energy.loc[charging_process.park_end-full_charging_steps:charging_process.park_end-1, cp] += \
-                charging_process.netto_charging_capacity
-            lower_energy.loc[charging_process.park_end - full_charging_steps-1, cp] += \
-                part_time_step * charging_process.netto_charging_capacity
+            lower_energy.loc[
+                charging_process.park_end
+                - full_charging_steps : charging_process.park_end
+                - 1,
+                cp,
+            ] += charging_process.netto_charging_capacity
+            lower_energy.loc[
+                charging_process.park_end - full_charging_steps - 1, cp
+            ] += (part_time_step * charging_process.netto_charging_capacity)
             # upper band
-            upper_energy.loc[charging_process.park_start-1: charging_process.park_start+full_charging_steps-2, cp] += \
-                charging_process.netto_charging_capacity
-            upper_energy.loc[charging_process.park_start+full_charging_steps-1, cp] += \
-                part_time_step * charging_process.netto_charging_capacity
+            upper_energy.loc[
+                charging_process.park_start
+                - 1 : charging_process.park_start
+                + full_charging_steps
+                - 2,
+                cp,
+            ] += charging_process.netto_charging_capacity
+            upper_energy.loc[
+                charging_process.park_start + full_charging_steps - 1, cp
+            ] += (part_time_step * charging_process.netto_charging_capacity)
     # sanity check
-    if ((lower_energy-upper_power * edisgo_obj.electromobility.eta_charging_points) > 1e-6).any().any():
-        raise ValueError("Lower energy has power values higher than nominal power. Please check.")
-    if ((upper_energy-upper_power * edisgo_obj.electromobility.eta_charging_points) > 1e-6).any().any():
-        raise ValueError("Upper energy has power values higher than nominal power. Please check.")
+    if (
+        (
+            (
+                lower_energy
+                - upper_power * edisgo_obj.electromobility.eta_charging_points
+            )
+            > 1e-6
+        )
+        .any()
+        .any()
+    ):
+        raise ValueError(
+            "Lower energy has power values higher than nominal power. Please check."
+        )
+    if (
+        (
+            (
+                upper_energy
+                - upper_power * edisgo_obj.electromobility.eta_charging_points
+            )
+            > 1e-6
+        )
+        .any()
+        .any()
+    ):
+        raise ValueError(
+            "Upper energy has power values higher than nominal power. Please check."
+        )
     if ((upper_energy.cumsum() - lower_energy.cumsum()) < -1e-6).any().any():
-        raise ValueError("Lower energy is higher than upper energy bound. Please check.")
+        raise ValueError(
+            "Lower energy is higher than upper energy bound. Please check."
+        )
     # Convert to MW and cumulate energy
     upper_power = upper_power / 1e3
     lower_energy = lower_energy.cumsum() / hourly_steps / 1e3
@@ -469,4 +505,8 @@ def _shorten_and_set_index(band):
     upper_power = _shorten_and_set_index(upper_power)
     lower_energy = _shorten_and_set_index(lower_energy)
     upper_energy = _shorten_and_set_index(upper_energy)
-    return upper_power, lower_energy, upper_energy
+    return {
+        "upper_power": upper_power,
+        "lower_energy": lower_energy,
+        "upper_energy": upper_energy,
+    }
diff --git a/edisgo/network/timeseries.py b/edisgo/network/timeseries.py
index 5b3664700..8660b2abe 100644
--- a/edisgo/network/timeseries.py
+++ b/edisgo/network/timeseries.py
@@ -1,16 +1,18 @@
-import logging
-import pandas as pd
-import numpy as np
 import datetime
+import logging
 import os
-
 from zipfile import ZipFile
+
+import numpy as np
+import pandas as pd
+from demandlib import bdew as bdew
+from demandlib import particular_profiles as profiles
 from workalendar.europe import Germany
-from demandlib import bdew as bdew, particular_profiles as profiles
-from edisgo.io.timeseries_import import import_feedin_timeseries
-from edisgo.tools.tools import assign_voltage_level_to_component,\
-    drop_duplicated_columns, get_weather_cells_intersecting_with_grid_district
 
+from edisgo.io.timeseries_import import import_feedin_timeseries
+from edisgo.tools.tools import (
+    assign_voltage_level_to_component, drop_duplicated_columns,
+    get_weather_cells_intersecting_with_grid_district)
 
 logger = logging.getLogger("edisgo")
 
@@ -29,10 +31,14 @@ def _get_attributes_to_save():
 
     """
     return [
-        "loads_active_power", "loads_reactive_power",
-        "generators_active_power", "generators_reactive_power",
-        "charging_points_active_power", "charging_points_reactive_power",
-        "storage_units_active_power", "storage_units_reactive_power"
+        "loads_active_power",
+        "loads_reactive_power",
+        "generators_active_power",
+        "generators_reactive_power",
+        "charging_points_active_power",
+        "charging_points_reactive_power",
+        "storage_units_active_power",
+        "storage_units_reactive_power",
     ]
 
 
@@ -288,8 +294,7 @@ def charging_points_reactive_power(self):
             return pd.DataFrame(index=self.timeindex)
 
     @charging_points_reactive_power.setter
-    def charging_points_reactive_power(self,
-                                       charging_points_reactive_power_ts):
+    def charging_points_reactive_power(self, charging_points_reactive_power_ts):
         self._charging_points_reactive_power = charging_points_reactive_power_ts
 
     # @property
@@ -359,10 +364,10 @@ def residual_load(self):
 
         """
         return (
-                self.loads_active_power.sum(axis=1) +
-                self.charging_points_active_power.sum(axis=1) -
-                self.generators_active_power.sum(axis=1) -
-                self.storage_units_active_power.sum(axis=1)
+            self.loads_active_power.sum(axis=1)
+            + self.charging_points_active_power.sum(axis=1)
+            - self.generators_active_power.sum(axis=1)
+            - self.storage_units_active_power.sum(axis=1)
         )
 
     @property
@@ -392,7 +397,7 @@ def timesteps_load_feedin_case(self):
         """
 
         return self.residual_load.apply(
-            lambda _: "feedin_case" if _ < 0. else "load_case"
+            lambda _: "feed-in_case" if _ < 0.0 else "load_case"
         )
 
     def reduce_memory(self, attr_to_reduce=None, to_type="float32"):
@@ -415,18 +420,17 @@ def reduce_memory(self, attr_to_reduce=None, to_type="float32"):
         """
         if attr_to_reduce is None:
             attr_to_reduce = [
-                "generators_active_power", "generators_reactive_power",
-                "loads_active_power", "loads_reactive_power",
+                "generators_active_power",
+                "generators_reactive_power",
+                "loads_active_power",
+                "loads_reactive_power",
                 "charging_points_active_power",
                 "charging_points_reactive_power",
-                "storage_units_active_power", "storage_units_reactive_power"
+                "storage_units_active_power",
+                "storage_units_reactive_power",
             ]
         for attr in attr_to_reduce:
-            setattr(
-                self,
-                attr,
-                getattr(self, attr).astype(to_type)
-            )
+            setattr(self, attr, getattr(self, attr).astype(to_type))
 
     def to_csv(self, directory, reduce_memory=False, **kwargs):
         """
@@ -470,8 +474,7 @@ def to_csv(self, directory, reduce_memory=False, **kwargs):
                     os.path.join(directory, "{}.csv".format(attr))
                 )
 
-    def from_csv(
-            self, data_path, dtype=None, from_zip_archive=False):
+    def from_csv(self, data_path, dtype=None, from_zip_archive=False):
         """
         Restores time series from csv files.
 
@@ -520,20 +523,15 @@ def from_csv(
             if from_zip_archive:
                 # open zip file to make it readable for pandas
                 with zip.open(file) as f:
-                    df = pd.read_csv(
-                        f, index_col=0, parse_dates=True, dtype=dtype)
+                    df = pd.read_csv(f, index_col=0, parse_dates=True, dtype=dtype)
             else:
                 path = os.path.join(data_path, file)
-                df = pd.read_csv(
-                    path, index_col=0, parse_dates=True, dtype=dtype)
+                df = pd.read_csv(path, index_col=0, parse_dates=True, dtype=dtype)
 
-            setattr(
-                self, attr, df)
+            setattr(self, attr, df)
 
             if timeindex is None:
-                timeindex = getattr(
-                    self,
-                    f"_{attr}").index
+                timeindex = getattr(self, f"_{attr}").index
 
         if from_zip_archive:
             # make sure to destroy ZipFile Class to close any open connections
@@ -639,8 +637,7 @@ def get_component_timeseries(edisgo_obj, **kwargs):
     mode = kwargs.get("mode", None)
     timeindex = kwargs.get("timeindex", edisgo_obj.timeseries.timeindex)
     # reset TimeSeries
-    edisgo_obj.timeseries = TimeSeries(
-        timeindex=timeindex)
+    edisgo_obj.timeseries = TimeSeries(timeindex=timeindex)
     edisgo_obj.timeseries.mode = mode
     if mode:
         if "worst-case" in mode:
@@ -656,48 +653,54 @@ def get_component_timeseries(edisgo_obj, **kwargs):
         elif mode == "manual":
             if kwargs.get("loads_active_power", None) is not None:
                 edisgo_obj.timeseries.loads_active_power = kwargs.get(
-                    "loads_active_power")
+                    "loads_active_power"
+                )
             if kwargs.get("loads_reactive_power", None) is not None:
                 edisgo_obj.timeseries.loads_reactive_power = kwargs.get(
-                    "loads_reactive_power")
+                    "loads_reactive_power"
+                )
 
             if kwargs.get("generators_active_power", None) is not None:
                 edisgo_obj.timeseries.generators_active_power = kwargs.get(
-                    "generators_active_power")
+                    "generators_active_power"
+                )
             if kwargs.get("generators_reactive_power", None) is not None:
                 edisgo_obj.timeseries.generators_reactive_power = kwargs.get(
-                    "generators_reactive_power")
+                    "generators_reactive_power"
+                )
 
             if kwargs.get("storage_units_active_power", None) is not None:
                 edisgo_obj.timeseries.storage_units_active_power = kwargs.get(
-                    "storage_units_active_power")
+                    "storage_units_active_power"
+                )
             if kwargs.get("storage_units_reactive_power", None) is not None:
-                edisgo_obj.timeseries.storage_units_reactive_power = \
-                    kwargs.get("storage_units_reactive_power")
+                edisgo_obj.timeseries.storage_units_reactive_power = kwargs.get(
+                    "storage_units_reactive_power"
+                )
 
             if kwargs.get("charging_points_active_power", None) is not None:
-                edisgo_obj.timeseries.charging_points_active_power = \
-                    kwargs.get("charging_points_active_power")
+                edisgo_obj.timeseries.charging_points_active_power = kwargs.get(
+                    "charging_points_active_power"
+                )
             if kwargs.get("charging_points_reactive_power", None) is not None:
-                edisgo_obj.timeseries.charging_points_reactive_power = \
-                    kwargs.get("charging_points_reactive_power")
+                edisgo_obj.timeseries.charging_points_reactive_power = kwargs.get(
+                    "charging_points_reactive_power"
+                )
         else:
             raise ValueError("{} is not a valid mode.".format(mode))
     else:
         config_data = edisgo_obj.config
 
-        weather_cell_ids = get_weather_cells_intersecting_with_grid_district(
-            edisgo_obj)
+        weather_cell_ids = get_weather_cells_intersecting_with_grid_district(edisgo_obj)
 
         # feed-in time series of fluctuating renewables
         ts = kwargs.get("timeseries_generation_fluctuating", None)
         if isinstance(ts, pd.DataFrame):
             edisgo_obj.timeseries.generation_fluctuating = ts
         elif isinstance(ts, str) and ts == "oedb":
-            edisgo_obj.timeseries.generation_fluctuating = \
-                import_feedin_timeseries(
-                    config_data, weather_cell_ids, kwargs.get(
-                        "timeindex", None))
+            edisgo_obj.timeseries.generation_fluctuating = import_feedin_timeseries(
+                config_data, weather_cell_ids, kwargs.get("timeindex", None)
+            )
         else:
             raise ValueError(
                 "Your input for "
@@ -739,8 +742,7 @@ def get_component_timeseries(edisgo_obj, **kwargs):
             )
         else:
             raise ValueError(
-                'Your input for "timeseries_load" is not '
-                "valid.".format(mode)
+                'Your input for "timeseries_load" is not ' "valid.".format(mode)
             )
         # reactive power timeseries for loads
         ts = kwargs.get("timeseries_load_reactive_power", None)
@@ -790,40 +792,38 @@ def _load_from_timeseries(edisgo_obj, load_names=None):
     )
     # set active power
     edisgo_obj.timeseries.loads_active_power = pd.concat(
-        [edisgo_obj.timeseries.loads_active_power,
-         loads.apply(
-             lambda x: edisgo_obj.timeseries.load[x.sector]
-             * x.annual_consumption
-             if x.sector in edisgo_obj.timeseries.load.columns
-             else edisgo_obj.timeseries.load['other']
-             * x.annual_consumption,
-             axis=1).T
-         ],
-        axis=1
+        [
+            edisgo_obj.timeseries.loads_active_power,
+            loads.apply(
+                lambda x: edisgo_obj.timeseries.load[x.sector] * x.annual_consumption
+                if x.sector in edisgo_obj.timeseries.load.columns
+                else edisgo_obj.timeseries.load["other"] * x.annual_consumption,
+                axis=1,
+            ).T,
+        ],
+        axis=1,
     )
 
     # if reactive power is given as attribute set with inserted timeseries
     if hasattr(edisgo_obj.timeseries, "load_reactive_power"):
         edisgo_obj.timeseries.loads_reactive_power = pd.concat(
-            [edisgo_obj.timeseries.loads_reactive_power,
-             loads.apply(
-                 lambda x: edisgo_obj.timeseries.load_reactive_power[x.sector]
-                 * x.annual_consumption
-                 if x.sector in
-                    edisgo_obj.timeseries.load_reactive_power.columns
-                 else edisgo_obj.timeseries.load_reactive_power['other']
-                 * x.annual_consumption,
-                 axis=1
-             )
-             ],
-            axis=1
+            [
+                edisgo_obj.timeseries.loads_reactive_power,
+                loads.apply(
+                    lambda x: edisgo_obj.timeseries.load_reactive_power[x.sector]
+                    * x.annual_consumption
+                    if x.sector in edisgo_obj.timeseries.load_reactive_power.columns
+                    else edisgo_obj.timeseries.load_reactive_power["other"]
+                    * x.annual_consumption,
+                    axis=1,
+                ),
+            ],
+            axis=1,
         )
     # set default reactive load
     else:
         _set_reactive_power_time_series_for_fixed_cosphi_using_config(
-            edisgo_obj=edisgo_obj,
-            df=loads,
-            component_type="loads"
+            edisgo_obj=edisgo_obj, df=loads, component_type="loads"
         )
 
 
@@ -841,20 +841,16 @@ def _timeseries_fluctuating():
             ).T
         else:
             return gens_fluctuating.apply(
-                lambda x: edisgo_obj.timeseries.generation_fluctuating[
-                    x.type
-                ].T
+                lambda x: edisgo_obj.timeseries.generation_fluctuating[x.type].T
                 * x.p_nom,
                 axis=1,
             ).T
 
     def _timeseries_dispatchable():
         return gens_dispatchable.apply(
-            lambda x: edisgo_obj.timeseries.generation_dispatchable[x.type]
-            * x.p_nom
+            lambda x: edisgo_obj.timeseries.generation_dispatchable[x.type] * x.p_nom
             if x.type in edisgo_obj.timeseries.generation_dispatchable.columns
-            else edisgo_obj.timeseries.generation_dispatchable["other"]
-            * x.p_nom,
+            else edisgo_obj.timeseries.generation_dispatchable["other"] * x.p_nom,
             axis=1,
         ).T
 
@@ -863,9 +859,7 @@ def _timeseries_dispatchable():
     # get all generators
     gens = edisgo_obj.topology.generators_df.loc[generator_names]
     # drop existing timeseries
-    _drop_existing_component_timeseries(
-        edisgo_obj, "generators", generator_names
-    )
+    _drop_existing_component_timeseries(edisgo_obj, "generators", generator_names)
     # handling of fluctuating generators
     gens_fluctuating = gens[gens.type.isin(["solar", "wind"])]
     gens_dispatchable = gens[~gens.index.isin(gens_fluctuating.index)]
@@ -879,7 +873,7 @@ def _timeseries_dispatchable():
                 _timeseries_dispatchable(),
             ],
             axis=1,
-            sort=False
+            sort=False,
         )
     if not gens_fluctuating.empty:
         edisgo_obj.timeseries.generators_active_power = pd.concat(
@@ -888,7 +882,7 @@ def _timeseries_dispatchable():
                 _timeseries_fluctuating(),
             ],
             axis=1,
-            sort=False
+            sort=False,
         )
 
     # set reactive power if given as attribute
@@ -900,18 +894,17 @@ def _timeseries_dispatchable():
     ):
 
         edisgo_obj.timeseries.generators_reactive_power = pd.concat(
-            [edisgo_obj.timeseries.generators_reactive_power,
-             edisgo_obj.timeseries.generation_reactive_power.loc[:, gens.index]
-             ],
-            axis=1
+            [
+                edisgo_obj.timeseries.generators_reactive_power,
+                edisgo_obj.timeseries.generation_reactive_power.loc[:, gens.index],
+            ],
+            axis=1,
         )
     # set default reactive power by cos_phi
     else:
         logger.debug("Reactive power calculated by cos(phi).")
         _set_reactive_power_time_series_for_fixed_cosphi_using_config(
-            edisgo_obj=edisgo_obj,
-            df=gens,
-            component_type="generators"
+            edisgo_obj=edisgo_obj, df=gens, component_type="generators"
         )
 
 
@@ -941,12 +934,8 @@ def _storage_from_timeseries(
     """
     if name_storage_units is None:
         name_storage_units = edisgo_obj.topology.storage_units_df.index
-    storage_units_df = edisgo_obj.topology.storage_units_df.loc[
-        name_storage_units
-    ]
-    _drop_existing_component_timeseries(
-        edisgo_obj, "storage_units", name_storage_units
-    )
+    storage_units_df = edisgo_obj.topology.storage_units_df.loc[name_storage_units]
+    _drop_existing_component_timeseries(edisgo_obj, "storage_units", name_storage_units)
 
     if len(storage_units_df) == 0:
         edisgo_obj.timeseries.storage_units_active_power = pd.DataFrame(
@@ -961,38 +950,39 @@ def _storage_from_timeseries(
     else:
         try:
             # check if indices and columns are correct
-            if (
-                ts_active_power.index == edisgo_obj.timeseries.timeindex
-            ).all():
-                edisgo_obj.timeseries.storage_units_active_power = drop_duplicated_columns(
-                    pd.concat(
-                        [edisgo_obj.timeseries.storage_units_active_power,
-                         ts_active_power.loc[:, name_storage_units]
-                         ],
-                        axis=1
+            if (ts_active_power.index == edisgo_obj.timeseries.timeindex).all():
+                edisgo_obj.timeseries.storage_units_active_power = (
+                    drop_duplicated_columns(
+                        pd.concat(
+                            [
+                                edisgo_obj.timeseries.storage_units_active_power,
+                                ts_active_power.loc[:, name_storage_units],
+                            ],
+                            axis=1,
+                        )
                     )
                 )
                 # check if reactive power is given
                 if (
                     ts_reactive_power is not None
-                    and (
-                        ts_active_power.index
-                        == edisgo_obj.timeseries.timeindex
-                    ).all()
+                    and (ts_active_power.index == edisgo_obj.timeseries.timeindex).all()
                 ):
-                    edisgo_obj.timeseries.storage_units_reactive_power = drop_duplicated_columns(
-                        pd.concat(
-                            [edisgo_obj.timeseries.storage_units_reactive_power,
-                             ts_reactive_power.loc[:, name_storage_units]
-                             ],
-                            axis=1
+                    edisgo_obj.timeseries.storage_units_reactive_power = (
+                        drop_duplicated_columns(
+                            pd.concat(
+                                [
+                                    edisgo_obj.timeseries.storage_units_reactive_power,
+                                    ts_reactive_power.loc[:, name_storage_units],
+                                ],
+                                axis=1,
+                            )
                         )
                     )
                 else:
                     _set_reactive_power_time_series_for_fixed_cosphi_using_config(
                         edisgo_obj=edisgo_obj,
                         df=storage_units_df,
-                        component_type="storage_units"
+                        component_type="storage_units",
                     )
             else:
                 raise ValueError(
@@ -1050,15 +1040,14 @@ def _worst_case_generation(edisgo_obj, modes, generator_names=None):
     worst_case_ts = pd.DataFrame(
         {
             "solar": [
-                worst_case_scale_factors["{}_feedin_pv".format(mode)]
-                for mode in modes
+                worst_case_scale_factors["{}_feed-in_pv".format(mode)] for mode in modes
             ],
             "wind": [
-                worst_case_scale_factors["{}_feedin_wind".format(mode)]
+                worst_case_scale_factors["{}_feed-in_wind".format(mode)]
                 for mode in modes
             ],
             "other": [
-                worst_case_scale_factors["{}_feedin_other".format(mode)]
+                worst_case_scale_factors["{}_feed-in_other".format(mode)]
                 for mode in modes
             ],
         },
@@ -1080,8 +1069,7 @@ def _worst_case_generation(edisgo_obj, modes, generator_names=None):
     cols_wind = gen_ts[gens_df.index[gens_df.type == "wind"]].columns
     if len(cols_wind) > 0:
         gen_ts[cols_wind] = pd.concat(
-            [worst_case_ts.loc[:, ["wind"]]] * len(cols_wind), axis=1,
-            sort=True
+            [worst_case_ts.loc[:, ["wind"]]] * len(cols_wind), axis=1, sort=True
         )
     # assign normalized active power time series to other generators
     cols = gen_ts.columns[~gen_ts.columns.isin(cols_pv.append(cols_wind))]
@@ -1091,22 +1079,17 @@ def _worst_case_generation(edisgo_obj, modes, generator_names=None):
         )
 
     # drop existing timeseries
-    _drop_existing_component_timeseries(
-        edisgo_obj, "generators", generator_names
-    )
+    _drop_existing_component_timeseries(edisgo_obj, "generators", generator_names)
 
     # multiply normalized time series by nominal power of generator
     edisgo_obj.timeseries.generators_active_power = pd.concat(
-        [edisgo_obj.timeseries.generators_active_power,
-         gen_ts.mul(gens_df.p_nom)],
-        axis=1
+        [edisgo_obj.timeseries.generators_active_power, gen_ts.mul(gens_df.p_nom)],
+        axis=1,
     )
 
     # calculate reactive power
     _set_reactive_power_time_series_for_fixed_cosphi_using_config(
-        edisgo_obj=edisgo_obj,
-        df=gens_df,
-        component_type="generators"
+        edisgo_obj=edisgo_obj, df=gens_df, component_type="generators"
     )
 
 
@@ -1140,18 +1123,14 @@ def _worst_case_load(edisgo_obj, modes, load_names=None):
     if check_loads.any():
         raise AttributeError(
             "The following loads have either missing bus, sector or "
-            "annual consumption: {}.".format(
-                check_loads[check_loads].index.values
-            )
+            "annual consumption: {}.".format(check_loads[check_loads].index.values)
         )
 
     # assign voltage level to loads
     if loads_df.empty:
         return
     loads_df["voltage_level"] = loads_df.apply(
-        lambda _: "lv"
-        if edisgo_obj.topology.buses_df.at[_.bus, "v_nom"] < 1
-        else "mv",
+        lambda _: "lv" if edisgo_obj.topology.buses_df.at[_.bus, "v_nom"] < 1 else "mv",
         axis=1,
     )
 
@@ -1171,10 +1150,7 @@ def _worst_case_load(edisgo_obj, modes, load_names=None):
     # assign power scaling factor to each load
     power_scaling_df = pd.DataFrame(
         data=np.transpose(
-            [
-                power_scaling[loads_df.at[col, "voltage_level"]]
-                for col in loads_df.index
-            ]
+            [power_scaling[loads_df.at[col, "voltage_level"]] for col in loads_df.index]
         ),
         index=edisgo_obj.timeseries.timeindex,
         columns=loads_df.index,
@@ -1187,16 +1163,15 @@ def _worst_case_load(edisgo_obj, modes, load_names=None):
 
     # calculate active power of loads
     edisgo_obj.timeseries.loads_active_power = pd.concat(
-        [edisgo_obj.timeseries.loads_active_power,
-         (power_scaling_df * loads_df.loc[:, "peak_load"])
-         ],
-        axis=1
+        [
+            edisgo_obj.timeseries.loads_active_power,
+            (power_scaling_df * loads_df.loc[:, "peak_load"]),
+        ],
+        axis=1,
     )
 
     _set_reactive_power_time_series_for_fixed_cosphi_using_config(
-        edisgo_obj=edisgo_obj,
-        df=loads_df,
-        component_type="loads"
+        edisgo_obj=edisgo_obj, df=loads_df, component_type="loads"
     )
 
 
@@ -1235,9 +1210,7 @@ def _worst_case_storage(edisgo_obj, modes, storage_names=None):
         if check_storage.any():
             raise AttributeError(
                 "The following storage units have either missing bus or "
-                "nominal power: {}.".format(
-                    check_storage[check_storage].index.values
-                )
+                "nominal power: {}.".format(check_storage[check_storage].index.values)
             )
 
         # active power
@@ -1260,18 +1233,17 @@ def _worst_case_storage(edisgo_obj, modes, storage_names=None):
         )
         edisgo_obj.timeseries.storage_units_active_power = drop_duplicated_columns(
             pd.concat(
-                [edisgo_obj.timeseries.storage_units_active_power,
-                (worst_case_ts * storage_df.p_nom)
-                 ],
-                axis=1
+                [
+                    edisgo_obj.timeseries.storage_units_active_power,
+                    (worst_case_ts * storage_df.p_nom),
+                ],
+                axis=1,
             ),
             keep="last",
         )
 
         _set_reactive_power_time_series_for_fixed_cosphi_using_config(
-            edisgo_obj=edisgo_obj,
-            df=storage_df,
-            component_type="storage_units"
+            edisgo_obj=edisgo_obj, df=storage_df, component_type="storage_units"
         )
 
 
@@ -1301,9 +1273,7 @@ def _check_timeindex(edisgo_obj):
             edisgo_obj.timeseries.storage_units_active_power.index
         ).all()
     except:
-        message = (
-            "Time index of feed-in and load time series does " "not match."
-        )
+        message = "Time index of feed-in and load time series does " "not match."
         logging.error(message)
         raise KeyError(message)
 
@@ -1343,9 +1313,7 @@ def add_loads_timeseries(edisgo_obj, load_names, **kwargs):
         if "worst-case" in edisgo_obj.timeseries.mode:
             modes = _get_worst_case_modes(edisgo_obj.timeseries.mode)
             # set random timeindex
-            _worst_case_load(
-                edisgo_obj=edisgo_obj, modes=modes, load_names=load_names
-            )
+            _worst_case_load(edisgo_obj=edisgo_obj, modes=modes, load_names=load_names)
         elif edisgo_obj.timeseries.mode == "manual":
             loads_active_power = kwargs.get("loads_active_power", None)
             if loads_active_power is not None:
@@ -1362,16 +1330,18 @@ def add_loads_timeseries(edisgo_obj, load_names, **kwargs):
             )
             # add new load timeseries
             edisgo_obj.timeseries.loads_active_power = pd.concat(
-                [edisgo_obj.timeseries.loads_active_power,
-                 loads_active_power.loc[:, load_names]
-                 ],
-                axis=1
+                [
+                    edisgo_obj.timeseries.loads_active_power,
+                    loads_active_power.loc[:, load_names],
+                ],
+                axis=1,
             )
             edisgo_obj.timeseries.loads_reactive_power = pd.concat(
-                [edisgo_obj.timeseries.loads_reactive_power,
-                 loads_reactive_power.loc[:, load_names]
-                 ],
-                axis=1
+                [
+                    edisgo_obj.timeseries.loads_reactive_power,
+                    loads_reactive_power.loc[:, load_names],
+                ],
+                axis=1,
             )
         else:
             raise ValueError(
@@ -1379,8 +1349,7 @@ def add_loads_timeseries(edisgo_obj, load_names, **kwargs):
             )
     else:
         # create load active and reactive power timeseries
-        _load_from_timeseries(
-            edisgo_obj=edisgo_obj, load_names=load_names)
+        _load_from_timeseries(edisgo_obj=edisgo_obj, load_names=load_names)
 
 
 def add_generators_timeseries(edisgo_obj, generator_names, **kwargs):
@@ -1446,33 +1415,34 @@ def add_generators_timeseries(edisgo_obj, generator_names, **kwargs):
             )
             # add new timeseries
             edisgo_obj.timeseries.generators_active_power = pd.concat(
-                [edisgo_obj.timeseries.generators_active_power,
-                 gens_active_power.loc[:, generator_names]
-                 ],
-                axis=1
+                [
+                    edisgo_obj.timeseries.generators_active_power,
+                    gens_active_power.loc[:, generator_names],
+                ],
+                axis=1,
             )
             edisgo_obj.timeseries.generators_reactive_power = pd.concat(
-                [edisgo_obj.timeseries.generators_reactive_power,
-                 gens_reactive_power.loc[:, generator_names]
-                 ],
-                axis=1
+                [
+                    edisgo_obj.timeseries.generators_reactive_power,
+                    gens_reactive_power.loc[:, generator_names],
+                ],
+                axis=1,
             )
         else:
             raise ValueError(
                 "{} is not a valid mode.".format(edisgo_obj.timeseries.mode)
             )
     else:
-        ts_dispatchable = kwargs.get(
-            "timeseries_generation_dispatchable", None
-        )
+        ts_dispatchable = kwargs.get("timeseries_generation_dispatchable", None)
         if ts_dispatchable is not None:
             if hasattr(edisgo_obj.timeseries, "generation_dispatchable"):
                 edisgo_obj.timeseries.generation_dispatchable = drop_duplicated_columns(
                     pd.concat(
-                        [edisgo_obj.timeseries.generation_dispatchable,
-                         ts_dispatchable
-                         ],
-                        axis=1
+                        [
+                            edisgo_obj.timeseries.generation_dispatchable,
+                            ts_dispatchable,
+                        ],
+                        axis=1,
                     ),
                     keep="last",
                 )
@@ -1482,19 +1452,20 @@ def add_generators_timeseries(edisgo_obj, generator_names, **kwargs):
         ts_reactive_power = kwargs.get("generation_reactive_power", None)
         if ts_reactive_power is not None:
             if hasattr(edisgo_obj.timeseries, "generation_reactive_power"):
-                edisgo_obj.timeseries.generation_reactive_power = drop_duplicated_columns(
-                    pd.concat(
-                        [edisgo_obj.timeseries.generation_reactive_power,
-                         ts_reactive_power
-                         ],
-                        axis=1
-                    ),
-                    keep="last",
-                )
-            else:
                 edisgo_obj.timeseries.generation_reactive_power = (
-                    ts_reactive_power
+                    drop_duplicated_columns(
+                        pd.concat(
+                            [
+                                edisgo_obj.timeseries.generation_reactive_power,
+                                ts_reactive_power,
+                            ],
+                            axis=1,
+                        ),
+                        keep="last",
+                    )
                 )
+            else:
+                edisgo_obj.timeseries.generation_reactive_power = ts_reactive_power
         # create load active and reactive power timeseries
         _generation_from_timeseries(
             edisgo_obj=edisgo_obj, generator_names=generator_names
@@ -1522,16 +1493,12 @@ def add_charging_points_timeseries(edisgo_obj, charging_point_names, **kwargs):
     """
     # TODO: only provision of time series is implemented, worst_case etc.
     #  is missing
-    ts_active_power = kwargs.get(
-        "ts_active_power", None
-    )
+    ts_active_power = kwargs.get("ts_active_power", None)
     if ts_active_power is not None:
         check_timeseries_for_index_and_cols(
             edisgo_obj, ts_active_power, charging_point_names
         )
-    ts_reactive_power = kwargs.get(
-        "ts_reactive_power", None
-    )
+    ts_reactive_power = kwargs.get("ts_reactive_power", None)
     if ts_reactive_power is not None:
         check_timeseries_for_index_and_cols(
             edisgo_obj,
@@ -1542,16 +1509,16 @@ def add_charging_points_timeseries(edisgo_obj, charging_point_names, **kwargs):
         edisgo_obj, "charging_points", charging_point_names
     )
     # add new timeseries
-    edisgo_obj.timeseries.charging_points_active_power = \
-        pd.concat([edisgo_obj.timeseries.charging_points_active_power,
-                   ts_active_power],
-                  axis=1, sort=False
-                  )
-    edisgo_obj.timeseries.charging_points_reactive_power = \
-        pd.concat([edisgo_obj.timeseries.charging_points_reactive_power,
-                   ts_reactive_power],
-                  axis=1, sort=False
-                  )
+    edisgo_obj.timeseries.charging_points_active_power = pd.concat(
+        [edisgo_obj.timeseries.charging_points_active_power, ts_active_power],
+        axis=1,
+        sort=False,
+    )
+    edisgo_obj.timeseries.charging_points_reactive_power = pd.concat(
+        [edisgo_obj.timeseries.charging_points_reactive_power, ts_reactive_power],
+        axis=1,
+        sort=False,
+    )
 
 
 def add_storage_units_timeseries(edisgo_obj, storage_unit_names, **kwargs):
@@ -1594,9 +1561,7 @@ def add_storage_units_timeseries(edisgo_obj, storage_unit_names, **kwargs):
                 storage_names=storage_unit_names,
             )
         elif edisgo_obj.timeseries.mode == "manual":
-            storage_units_active_power = kwargs.get(
-                "storage_units_active_power", None
-            )
+            storage_units_active_power = kwargs.get("storage_units_active_power", None)
             if storage_units_active_power is not None:
                 check_timeseries_for_index_and_cols(
                     edisgo_obj, storage_units_active_power, storage_unit_names
@@ -1615,16 +1580,18 @@ def add_storage_units_timeseries(edisgo_obj, storage_unit_names, **kwargs):
             )
             # add new storage timeseries
             edisgo_obj.timeseries.storage_units_active_power = pd.concat(
-                [edisgo_obj.timeseries.storage_units_active_power,
-                 storage_units_active_power.loc[:, storage_unit_names]
-                 ],
-                axis=1
+                [
+                    edisgo_obj.timeseries.storage_units_active_power,
+                    storage_units_active_power.loc[:, storage_unit_names],
+                ],
+                axis=1,
             )
             edisgo_obj.timeseries.storage_units_reactive_power = pd.concat(
-                [edisgo_obj.timeseries.storage_units_reactive_power,
-                 storage_units_reactive_power.loc[:, storage_unit_names]
-                 ],
-                axis=1
+                [
+                    edisgo_obj.timeseries.storage_units_reactive_power,
+                    storage_units_reactive_power.loc[:, storage_unit_names],
+                ],
+                axis=1,
             )
         else:
             raise ValueError(
@@ -1665,9 +1632,7 @@ def _drop_existing_component_timeseries(edisgo_obj, comp_type, comp_names):
         edisgo_obj.timeseries,
         comp_type + "_active_power",
         getattr(edisgo_obj.timeseries, comp_type + "_active_power").drop(
-            getattr(
-                edisgo_obj.timeseries, comp_type + "_active_power"
-            ).columns[
+            getattr(edisgo_obj.timeseries, comp_type + "_active_power").columns[
                 getattr(
                     edisgo_obj.timeseries, comp_type + "_active_power"
                 ).columns.isin(comp_names)
@@ -1679,9 +1644,7 @@ def _drop_existing_component_timeseries(edisgo_obj, comp_type, comp_names):
         edisgo_obj.timeseries,
         comp_type + "_reactive_power",
         getattr(edisgo_obj.timeseries, comp_type + "_reactive_power").drop(
-            getattr(
-                edisgo_obj.timeseries, comp_type + "_reactive_power"
-            ).columns[
+            getattr(edisgo_obj.timeseries, comp_type + "_reactive_power").columns[
                 getattr(
                     edisgo_obj.timeseries, comp_type + "_reactive_power"
                 ).columns.isin(comp_names)
@@ -1691,9 +1654,7 @@ def _drop_existing_component_timeseries(edisgo_obj, comp_type, comp_names):
     )
 
 
-def check_timeseries_for_index_and_cols(
-    edisgo_obj, timeseries, component_names
-):
+def check_timeseries_for_index_and_cols(edisgo_obj, timeseries, component_names):
     """
     Checks index and column names of inserted timeseries to make sure, they
     have the right format.
@@ -1783,9 +1744,7 @@ def _load_timeseries_demandlib(config_data, year):
         elec_demand = e_slp.get_profile(sectoral_consumption)
 
         # Add the slp for the industrial group
-        ilp = profiles.IndustrialLoadProfile(
-            e_slp.date_time_index, holidays=holidays
-        )
+        ilp = profiles.IndustrialLoadProfile(e_slp.date_time_index, holidays=holidays)
 
         # Beginning and end of workday, weekdays and weekend days, and scaling
         # factors by default
@@ -1825,9 +1784,7 @@ def _load_timeseries_demandlib(config_data, year):
                 raise TypeError
         except TypeError:
             year = datetime.date.today().year - 1
-            logger.warning(
-                "No valid year inserted. Year set to previous year."
-            )
+            logger.warning("No valid year inserted. Year set to previous year.")
         load = _load_timeseries_demandlib(config_data, year)
         load.rename(
             columns={
@@ -1858,13 +1815,14 @@ def _get_worst_case_modes(mode):
         list which can contains 'feedin-case', 'load_case' or both
     """
     if mode == "worst-case":
-        modes = ["feedin_case", "load_case"]
-    elif mode == "worst-case-feedin" or mode == "worst-case-load":
+        modes = ["feed-in_case", "load_case"]
+    elif mode == "worst-case-feed-in" or mode == "worst-case-load":
         modes = ["{}_case".format(mode.split("-")[-1])]
     else:
         raise ValueError("{} is not a valid mode.".format(mode))
     return modes
 
+
 def _get_q_sign_generator(reactive_power_mode):
     """
     Get the sign of reactive power in generator sign convention.
@@ -1956,7 +1914,8 @@ def fixed_cosphi(active_power, q_sign, power_factor):
 
 
 def _set_reactive_power_time_series_for_fixed_cosphi_using_config(
-        edisgo_obj, df, component_type):
+    edisgo_obj, df, component_type
+):
     """
     Calculates reactive power in Mvar for a fixed cosphi operation.
 
@@ -2013,37 +1972,34 @@ def _set_reactive_power_time_series_for_fixed_cosphi_using_config(
     else:
         raise ValueError(
             "Given 'component_type' is not valid. Valid options are "
-            "'generators','storage_units' and 'loads'.")
+            "'generators','storage_units' and 'loads'."
+        )
     for voltage_level in voltage_levels:
         cols = df.index[df.voltage_level == voltage_level]
         if len(cols) > 0:
             q_sign[cols] = get_q_sign(
-                reactive_power_mode[
-                    "{}_gen".format(voltage_level)
-                ]
+                reactive_power_mode["{}_gen".format(voltage_level)]
             )
-            power_factor[cols] = reactive_power_factor[
-                "{}_gen".format(voltage_level)
-            ]
+            power_factor[cols] = reactive_power_factor["{}_gen".format(voltage_level)]
 
     # calculate reactive power time series and append to TimeSeries object
     reactive_power_df = drop_duplicated_columns(
         pd.concat(
-            [getattr(edisgo_obj.timeseries,
-                     component_type + "_reactive_power"),
-             fixed_cosphi(
-                 getattr(edisgo_obj.timeseries,
-                         component_type + "_active_power").loc[:, df.index],
-                 q_sign,
-                 power_factor
-             )],
-            axis=1
+            [
+                getattr(edisgo_obj.timeseries, component_type + "_reactive_power"),
+                fixed_cosphi(
+                    getattr(
+                        edisgo_obj.timeseries, component_type + "_active_power"
+                    ).loc[:, df.index],
+                    q_sign,
+                    power_factor,
+                ),
+            ],
+            axis=1,
         ),
-        keep="last"
+        keep="last",
     )
 
     setattr(
-        edisgo_obj.timeseries,
-        component_type + "_reactive_power",
-        reactive_power_df
+        edisgo_obj.timeseries, component_type + "_reactive_power", reactive_power_df
     )
diff --git a/edisgo/opf/lopf.py b/edisgo/opf/lopf.py
index a0e118b51..bfa8af7cd 100644
--- a/edisgo/opf/lopf.py
+++ b/edisgo/opf/lopf.py
@@ -1,5 +1,3 @@
-# fmt: off
-
 # Methods to perform linearised DistFlow
 import itertools
 from copy import deepcopy
@@ -49,7 +47,7 @@ def prepare_time_invariant_parameters(
     pu=True,
     optimize_storage=True,
     optimize_ev_charging=True,
-    **kwargs
+    **kwargs,
 ):
     """
     Prepare parameters that do not change within the iterations of the rolling horizon
@@ -187,229 +185,13 @@ def prepare_time_invariant_parameters(
     return parameters
 
 
-def setup_grid_object(object):
-    """
-    Set up the grid and edisgo object.
-    """
-    if hasattr(object, "topology"):
-        grid_object = deepcopy(object.topology)
-        edisgo_object = deepcopy(object)
-        slack = grid_object.mv_grid.station.index
-    else:
-        grid_object = deepcopy(object)
-        edisgo_object = deepcopy(object.edisgo_obj)
-        slack = [
-            grid_object.transformers_df.bus1.iloc[0]
-        ]  # Todo: careful with MV grid, does not work with that right?
-    return edisgo_object, grid_object, slack
-
-
-def concat_parallel_branch_elements(grid_object):
-    """
-    Method to merge parallel lines and transformers into one element, respectively.
-
-    Parameters
-    ----------
-    grid_object
-
-    Returns
-    -------
-
-    """
-    lines = fuse_parallel_branches(grid_object.lines_df)
-    trafos = grid_object.transformers_df.loc[
-        grid_object.transformers_df.bus0.isin(grid_object.buses_df.index)
-    ].loc[grid_object.transformers_df.bus1.isin(grid_object.buses_df.index)]
-    transformers = fuse_parallel_branches(trafos)
-    return pd.concat([lines, transformers], sort=False)
-
-
-def fuse_parallel_branches(branches):
-    branches_tmp = branches[["bus0", "bus1"]]
-    parallel_branches = pd.DataFrame(columns=branches.columns)
-    if branches_tmp.duplicated().any():
-        duplicated_branches = branches_tmp.loc[branches_tmp.duplicated(keep=False)]
-        duplicated_branches["visited"] = False
-        branches_tmp.drop(duplicated_branches.index, inplace=True)
-        for name, buses in duplicated_branches.iterrows():
-            if duplicated_branches.loc[name, "visited"]:
-                continue
-            else:
-                parallel_branches_tmp = duplicated_branches.loc[
-                    (duplicated_branches == buses).all(axis=1)
-                ]
-                duplicated_branches.loc[parallel_branches_tmp.index, "visited"] = True
-                name_par = "_".join(str.split(name, "_")[:-1])
-                parallel_branches.loc[name_par] = branches.loc[name]
-                parallel_branches.loc[
-                    name_par, ["r", "x", "s_nom"]
-                ] = calculate_impedance_for_parallel_components(
-                    branches.loc[parallel_branches_tmp.index, ["r", "x", "s_nom"]],
-                    pu=False,
-                )
-    fused_branches = pd.concat(
-        [branches.loc[branches_tmp.index], parallel_branches], sort=False
-    )
-    return fused_branches
-
-
-def get_underlying_elements(parameters):
-    def _get_underlying_elements(
-        downstream_elements, power_factors, parameters, branch
-    ):
-        bus0 = parameters["branches"].loc[branch, "bus0"]
-        bus1 = parameters["branches"].loc[branch, "bus1"]
-        s_nom = parameters["branches"].loc[branch, parameters["pars"]["s_nom"]]
-        relevant_buses_bus0 = (
-            parameters["downstream_nodes_matrix"]
-            .loc[bus0][parameters["downstream_nodes_matrix"].loc[bus0] == 1]
-            .index.values
-        )
-        relevant_buses_bus1 = (
-            parameters["downstream_nodes_matrix"]
-            .loc[bus1][parameters["downstream_nodes_matrix"].loc[bus1] == 1]
-            .index.values
-        )
-        relevant_buses = list(
-            set(relevant_buses_bus0).intersection(relevant_buses_bus1)
-        )
-        downstream_elements.loc[branch, "buses"] = relevant_buses
-        if (
-            parameters["nodal_reactive_power"]
-            .loc[relevant_buses]
-            .sum()
-            .divide(s_nom)
-            .apply(abs)
-            > 1
-        ).any():
-            print(
-                "Careful: Reactive power already exceeding line capacity for branch "
-                "{}.".format(branch)
-            )
-        power_factors.loc[branch] = (
-            1
-            - parameters["nodal_reactive_power"]
-            .loc[relevant_buses]
-            .sum()
-            .divide(s_nom)
-            .apply(np.square)
-        ).apply(np.sqrt)
-        if "optimized_storage_units" in parameters:
-            downstream_elements.loc[branch, "flexible_storage"] = (
-                parameters["grid_object"]
-                .storage_units_df.loc[
-                    parameters["grid_object"].storage_units_df.index.isin(
-                        parameters["optimized_storage_units"]
-                    )
-                    & parameters["grid_object"].storage_units_df.bus.isin(
-                        relevant_buses
-                    )
-                ]
-                .index.values
-            )
-        else:
-            downstream_elements.loc[branch, "flexible_storage"] = []
-        if "optimized_charging_points" in parameters:
-            downstream_elements.loc[branch, "flexible_ev"] = (
-                parameters["grid_object"]
-                .charging_points_df.loc[
-                    parameters["grid_object"].charging_points_df.index.isin(
-                        parameters["optimized_charging_points"]
-                    )
-                    & parameters["grid_object"].charging_points_df.bus.isin(
-                        relevant_buses
-                    )
-                ]
-                .index.values
-            )
-        else:
-            downstream_elements.loc[branch, "flexible_ev"] = []
-        return downstream_elements, power_factors
-
-    downstream_elements = pd.DataFrame(
-        index=parameters["branches"].index,
-        columns=["buses", "flexible_storage", "flexible_ev"],
-    )
-    power_factors = pd.DataFrame(
-        index=parameters["branches"].index,
-        columns=parameters["nodal_active_power"].columns,
-    )
-    for branch in downstream_elements.index:
-        downstream_elements, power_factors = _get_underlying_elements(
-            downstream_elements, power_factors, parameters, branch
-        )
-    if power_factors.isna().any().any():
-        print(
-            "WARNING: Branch {} is overloaded with reactive power. Still needs "
-            "handling.".format(branch)
-        )
-        power_factors = power_factors.fillna(
-            0.01
-        )  # Todo: ask Gaby and Birgit about this
-    return downstream_elements, power_factors
-
-
-def get_residual_load_of_not_optimized_components(
-    grid,
-    edisgo,
-    relevant_storage_units=None,
-    relevant_charging_points=None,
-    relevant_generators=None,
-    relevant_loads=None,
-):
-    """
-    Method to get residual load of fixed components.
-
-    Parameters
-    ----------
-    grid
-    edisgo
-    relevant_storage_units
-    relevant_charging_points
-    relevant_generators
-    relevant_loads
-
-    Returns
-    -------
-
-    """
-    if relevant_loads is None:
-        relevant_loads = grid.loads_df.index
-    if relevant_generators is None:
-        relevant_generators = grid.generators_df.index
-    if relevant_storage_units is None:
-        relevant_storage_units = grid.storage_units_df.index
-    if relevant_charging_points is None:
-        relevant_charging_points = grid.charging_points_df.index
-
-    if edisgo.timeseries.charging_points_active_power.empty:
-        return (
-            edisgo.timeseries.generators_active_power[relevant_generators].sum(axis=1)
-            + edisgo.timeseries.storage_units_active_power[relevant_storage_units].sum(
-                axis=1
-            )
-            - edisgo.timeseries.loads_active_power[relevant_loads].sum(axis=1)
-        ).loc[edisgo.timeseries.timeindex]
-    else:
-        return (
-            edisgo.timeseries.generators_active_power[relevant_generators].sum(axis=1)
-            + edisgo.timeseries.storage_units_active_power[relevant_storage_units].sum(
-                axis=1
-            )
-            - edisgo.timeseries.loads_active_power[relevant_loads].sum(axis=1)
-            - edisgo.timeseries.charging_points_active_power[
-                relevant_charging_points
-            ].sum(axis=1)
-        ).loc[edisgo.timeseries.timeindex]
-
-
 def setup_model(
     timeinvariant_parameters,
     timesteps,
     optimize_storage=True,
     optimize_ev_charging=True,
     objective="curtailment",
-    **kwargs
+    **kwargs,
 ):
     """
     Method to set up pyomo model for optimisation of storage procurement
@@ -426,84 +208,28 @@ def setup_model(
     :return:
     """
 
-    def init_active_nodal_power(model, bus, time):
-        return (
-            timeinvariant_parameters["nodal_active_power"]
-            .T.loc[model.timeindex[time]]
-            .loc[bus]
-        )
-
-    def init_reactive_nodal_power(model, bus, time):
-        return (
-            timeinvariant_parameters["nodal_reactive_power"]
-            .T.loc[model.timeindex[time]]
-            .loc[bus]
-        )
-
-    def init_active_nodal_load(model, bus, time):
-        return (
-            timeinvariant_parameters["nodal_active_load"]
-            .T.loc[model.timeindex[time]]
-            .loc[bus]
-        )
-
-    def init_reactive_nodal_load(model, bus, time):
-        return (
-            timeinvariant_parameters["nodal_reactive_load"]
-            .T.loc[model.timeindex[time]]
-            .loc[bus]
-        )
-
-    def init_active_nodal_feedin(model, bus, time):
-        return (
-            timeinvariant_parameters["nodal_active_feedin"]
-            .T.loc[model.timeindex[time]]
-            .loc[bus]
-        )
-
-    def init_reactive_nodal_feedin(model, bus, time):
-        return (
-            timeinvariant_parameters["nodal_reactive_feedin"]
-            .T.loc[model.timeindex[time]]
-            .loc[bus]
-        )
-
-    def init_power_factors(model, branch, time):
-        return timeinvariant_parameters["power_factors"].loc[
-            branch, model.timeindex[time]
-        ]
-
-    t1 = perf_counter()
-    model = pm.ConcreteModel()
-    # check if correct value of objective is inserted
-    if objective not in [
-        "curtailment",
-        "peak_load",
-        "minimize_energy_level",
-        "residual_load",
-        "maximize_energy_level",
-        "minimize_loading",
-    ]:
-        raise ValueError("The objective you inserted is not implemented yet.")
-
-    edisgo_object, grid_object, slack = (
-        timeinvariant_parameters["edisgo_object"],
-        timeinvariant_parameters["grid_object"],
-        timeinvariant_parameters["slack"],
-    )
-    # check if multiple voltage levels are present
-    if len(grid_object.buses_df.v_nom.unique()) > 1:
-        print(
-            "More than one voltage level included. Please make sure to "
-            "adapt all impedance values to one reference system."
-        )
-
-    # Todo: Extract kwargs values from cfg?
+    # Todo: Extract kwargs values from cfg?
+    t1 = perf_counter()
+    model = pm.ConcreteModel()
+    # check if correct value of objective is inserted
+    if objective not in [
+        "curtailment",
+        "peak_load",
+        "minimize_energy_level",
+        "residual_load",
+        "maximize_energy_level",
+        "minimize_loading",
+    ]:
+        raise ValueError("The objective you inserted is not implemented yet.")
+    # Todo: unnecessary?
+    edisgo_object, grid_object, slack = (
+        timeinvariant_parameters["edisgo_object"],
+        timeinvariant_parameters["grid_object"],
+        timeinvariant_parameters["slack"],
+    )
 
     # DEFINE SETS AND FIX PARAMETERS
     print("Setup model: Defining sets and parameters.")
-    model.bus_set = pm.Set(initialize=grid_object.buses_df.index)
-    model.slack_bus = pm.Set(initialize=slack)
     model.time_set = pm.RangeSet(0, len(timesteps) - 1)
     model.time_zero = [model.time_set.at(1)]
     overlap_interations = kwargs.get("overlap_interations", None)
@@ -511,6 +237,7 @@ def init_power_factors(model, branch, time):
         model.time_end = [model.time_set.at(-overlap_interations)]
     else:
         model.time_final = [model.time_set.at(-1)]
+        model.time_end = [model.time_set.at(-1)]
     model.time_non_zero = model.time_set - [model.time_set.at(1)]
     model.times_fixed_soc = pm.Set(
         initialize=[model.time_set.at(1), model.time_set.at(-1)]
@@ -532,90 +259,12 @@ def init_power_factors(model, branch, time):
         )
         model.fixed_storage_set = model.storage_set - model.optimized_storage_set
         model.fix_relative_soc = kwargs.get("fix_relative_soc", 0.5)
-    if optimize_ev_charging:
-        model.charging_points_set = pm.Set(
-            initialize=grid_object.charging_points_df.index
-        )
-        model.flexible_charging_points_set = pm.Set(
-            initialize=timeinvariant_parameters["optimized_charging_points"]
-        )
-        model.inflexible_charging_points_set = (
-            model.charging_points_set - model.flexible_charging_points_set
-        )
-        model.upper_ev_power = timeinvariant_parameters["ev_flex_bands"]["upper_power"]
-        model.upper_ev_energy = timeinvariant_parameters["ev_flex_bands"][
-            "upper_energy"
-        ]
-        model.lower_ev_energy = timeinvariant_parameters["ev_flex_bands"][
-            "lower_energy"
-        ]
-        model.charging_efficiency = kwargs.get("charging_efficiency", 0.9)
-        model.lower_bound_ev = pm.Param(
-            model.flexible_charging_points_set,
-            model.time_set,
-            initialize=set_lower_band_ev,
-            mutable=True,
-        )
-        model.upper_bound_ev = pm.Param(
-            model.flexible_charging_points_set,
-            model.time_set,
-            initialize=set_upper_band_ev,
-            mutable=True,
-        )
-        model.power_bound_ev = pm.Param(
-            model.flexible_charging_points_set,
-            model.time_set,
-            initialize=set_power_band_ev,
-            mutable=True,
-        )
-    model.v_min = kwargs.get("v_min", 0.9)
-    model.v_max = kwargs.get("v_max", 1.1)
-    model.v_nom = timeinvariant_parameters["v_nom"]
-    model.thermal_limit = kwargs.get("thermal_limit", 1.0)
-    model.pars = timeinvariant_parameters["pars"]
+
     res_load = {
         i: timeinvariant_parameters["res_load_inflexible_units"][model.timeindex[i]]
         for i in model.time_set
     }
     model.residual_load = pm.Param(model.time_set, initialize=res_load, mutable=True)
-    model.grid = grid_object
-    model.downstream_nodes_matrix = timeinvariant_parameters["downstream_nodes_matrix"]
-
-    model.nodal_active_power = pm.Param(
-        model.bus_set, model.time_set, initialize=init_active_nodal_power, mutable=True
-    )
-    model.nodal_reactive_power = pm.Param(
-        model.bus_set,
-        model.time_set,
-        initialize=init_reactive_nodal_power,
-        mutable=True,
-    )
-    model.nodal_active_load = pm.Param(
-        model.bus_set, model.time_set, initialize=init_active_nodal_load, mutable=True
-    )
-    model.nodal_reactive_load = pm.Param(
-        model.bus_set, model.time_set, initialize=init_reactive_nodal_load, mutable=True
-    )
-    model.nodal_active_feedin = pm.Param(
-        model.bus_set, model.time_set, initialize=init_active_nodal_feedin, mutable=True
-    )
-    model.nodal_reactive_feedin = pm.Param(
-        model.bus_set,
-        model.time_set,
-        initialize=init_reactive_nodal_feedin,
-        mutable=True,
-    )
-    model.tan_phi_load = timeinvariant_parameters["tan_phi_load"]
-    model.tan_phi_feedin = timeinvariant_parameters["tan_phi_feedin"]
-    model.v_slack = kwargs.get("v_slack", model.v_nom)
-    model.branches = timeinvariant_parameters["branches"]
-    model.branch_set = pm.Set(initialize=model.branches.index)
-    model.underlying_branch_elements = timeinvariant_parameters[
-        "underlying_branch_elements"
-    ]
-    model.power_factors = pm.Param(
-        model.branch_set, model.time_set, initialize=init_power_factors, mutable=True
-    )
 
     if objective == "peak_load":
         model.delta_min = kwargs.get("delta_min", 0.9)
@@ -625,53 +274,8 @@ def init_power_factors(model, branch, time):
     elif objective == "minimize_energy_level" or objective == "maximize_energy_level":
         model.grid_power_flexible = pm.Var(model.time_set)
 
-    # add n-1 security
-    # adapt i_lines_allowed for radial feeders
-    buses_in_cycles = list(
-        set(itertools.chain.from_iterable(edisgo_object.topology.rings))
-    )
-
-    # Find lines in cycles
-    lines_in_cycles = list(
-        grid_object.lines_df.loc[
-            grid_object.lines_df[["bus0", "bus1"]].isin(buses_in_cycles).all(axis=1)
-        ].index.values
-    )
-
-    model.branches_load_factors = pd.DataFrame(
-        index=model.time_set, columns=model.branch_set
-    )
-    model.branches_load_factors.loc[:, :] = 1
-    tmp_residual_load = edisgo_object.timeseries.residual_load.loc[timesteps]
-    indices = pd.DataFrame(index=timesteps, columns=["index"])
-    indices["index"] = [i for i in range(len(timesteps))]
-    model.branches_load_factors.loc[
-        indices.loc[tmp_residual_load.loc[timesteps] < 0].values.T[0], lines_in_cycles
-    ] = kwargs.get(
-        "load_factor_rings", 1.0
-    )  # 0.5
-
     # DEFINE VARIABLES
-    print("Setup model: Defining variables.")
-    model.p_cum = pm.Var(model.branch_set, model.time_set)
-    model.slack_p_cum_pos = pm.Var(model.branch_set, model.time_set, bounds=(0, None))
-    model.slack_p_cum_neg = pm.Var(model.branch_set, model.time_set, bounds=(0, None))
-    model.q_cum = pm.Var(model.branch_set, model.time_set)
-    model.v = pm.Var(model.bus_set, model.time_set)
-    model.slack_v_pos = pm.Var(model.bus_set, model.time_set, bounds=(0, None))
-    model.slack_v_neg = pm.Var(model.bus_set, model.time_set, bounds=(0, None))
-    # if not objective == 'minimize_energy_level' and \
-    #         not objective == 'maximize_energy_level':
-    model.curtailment_load = pm.Var(
-        model.bus_set,
-        model.time_set,
-        bounds=lambda m, b, t: (0, m.nodal_active_load[b, t]),
-    )
-    model.curtailment_feedin = pm.Var(
-        model.bus_set,
-        model.time_set,
-        bounds=lambda m, b, t: (0, m.nodal_active_feedin[b, t]),
-    )
+
     if optimize_storage:
         model.soc = pm.Var(
             model.optimized_storage_set,
@@ -689,51 +293,40 @@ def init_power_factors(model, branch, time):
                 m.grid.storage_units_df.loc[b, model.pars["p_nom"]],
             ),
         )
-    if optimize_ev_charging:
 
-        model.charging_ev = pm.Var(
-            model.flexible_charging_points_set,
-            model.time_set,
-            bounds=lambda m, b, t: (0, m.power_bound_ev[b, t]),
+    if optimize_ev_charging:
+        print("Setup model: Adding EV model.")
+        model = add_ev_model_bands(
+            model=model,
+            timeinvariant_parameters=timeinvariant_parameters,
+            grid_object=grid_object,
+            charging_efficiency=kwargs.get("charging_efficiency", 0.9),
+            energy_level_start=kwargs.get("energy_level_start", None),
+            energy_level_end=kwargs.get("energy_level_end", None),
+            energy_level_beginning=kwargs.get("energy_level_beginning", None),
+            charging_start=kwargs.get("charging_start", None),
+        )
+
+    if not objective == "minimize_energy_level" or objective == "maximize_energy_level":
+        print("Setup model: Adding grid model.")
+        model = add_grid_model_lopf(
+            model=model,
+            timeinvariant_parameters=timeinvariant_parameters,
+            timesteps=timesteps,
+            edisgo_object=edisgo_object,
+            grid_object=grid_object,
+            slack=slack,
+            v_min=kwargs.get("v_min", 0.9),
+            v_max=kwargs.get("v_max", 1.1),
+            thermal_limits=kwargs.get("thermal_limit", 1.0),
+            v_slack=kwargs.get("v_slack", timeinvariant_parameters["v_nom"]),
+            load_factor_rings=kwargs.get("load_factor_rings", 1.0),
+            # 0.5 Todo: change to edisgo.config["grid_expansion_load_factors"]
+            #  ["mv_load_case_line"]?
         )
 
-        model.curtailment_ev = pm.Var(model.bus_set, model.time_set, bounds=(0, None))
-        if not (
-            objective == "minimize_energy_level" or objective == "maximize_energy_level"
-        ):
-            model.energy_level_ev = pm.Var(
-                model.flexible_charging_points_set,
-                model.time_set,
-                bounds=lambda m, b, t: (m.lower_bound_ev[b, t], m.upper_bound_ev[b, t]),
-            )
-
     # DEFINE CONSTRAINTS
-    print("Setup model: Setting constraints.")
-    model.ActivePower = pm.Constraint(
-        model.branch_set, model.time_set, rule=active_power
-    )
-    model.UpperActive = pm.Constraint(
-        model.branch_set, model.time_set, rule=upper_active_power
-    )
-    model.LowerActive = pm.Constraint(
-        model.branch_set, model.time_set, rule=lower_active_power
-    )
-    # model.ReactivePower = pm.Constraint(model.branch_set, model.time_set,
-    #                                     rule=reactive_power)
-    model.SlackVoltage = pm.Constraint(
-        model.slack_bus, model.time_set, rule=slack_voltage
-    )
-    model.VoltageDrop = pm.Constraint(
-        model.branch_set, model.time_set, rule=voltage_drop
-    )
-    model.UpperVoltage = pm.Constraint(
-        model.bus_set, model.time_set, rule=upper_voltage
-    )
-    model.LowerVoltage = pm.Constraint(
-        model.bus_set, model.time_set, rule=lower_voltage
-    )
-    # model.UpperCurtLoad = pm.Constraint(model.bus_set, model.time_set,
-    #                                     rule=upper_bound_curtailment_load)
+
     if optimize_storage:
         model.BatteryCharging = pm.Constraint(
             model.storage_set, model.time_non_zero, rule=soc
@@ -741,99 +334,6 @@ def init_power_factors(model, branch, time):
         model.FixedSOC = pm.Constraint(
             model.storage_set, model.times_fixed_soc, rule=fix_soc
         )
-    if optimize_ev_charging and not (
-        objective == "minimize_energy_level" or objective == "maximize_energy_level"
-    ):
-        model.EVCharging = pm.Constraint(
-            model.flexible_charging_points_set, model.time_non_zero, rule=charging_ev
-        )
-        model.UpperCurtEV = pm.Constraint(
-            model.bus_set, model.time_set, rule=upper_bound_curtailment_ev
-        )
-        # set initial energy level
-        energy_level_start = kwargs.get("energy_level_start", None)
-        model.energy_level_start = pm.Param(
-            model.flexible_charging_points_set,
-            initialize=energy_level_start,
-            mutable=True,
-        )
-        model.slack_initial_energy_pos = pm.Var(
-            model.flexible_charging_points_set, bounds=(0, None)
-        )
-        model.slack_initial_energy_neg = pm.Var(
-            model.flexible_charging_points_set, bounds=(0, None)
-        )
-        model.InitialEVEnergyLevel = pm.Constraint(
-            model.flexible_charging_points_set,
-            model.time_zero,
-            rule=initial_energy_level,
-        )
-        model.InitialEVEnergyLevelStart = pm.Constraint(
-            model.flexible_charging_points_set, model.time_zero, rule=fixed_energy_level
-        )
-        if energy_level_start is None:
-            model.InitialEVEnergyLevel.deactivate()
-        else:
-            model.InitialEVEnergyLevelStart.deactivate()
-        # set final energy level and if necessary charging power
-        energy_level_end = kwargs.get("energy_level_end", None)
-        model.energy_level_end = pm.Param(
-            model.flexible_charging_points_set,
-            initialize=energy_level_end,
-            mutable=True,
-        )
-        model.FinalEVEnergyLevelFix = pm.Constraint(
-            model.flexible_charging_points_set, model.time_end, rule=fixed_energy_level
-        )
-
-        energy_level_beginning = kwargs.get("energy_level_beginning", None)
-        if energy_level_beginning is None:
-            model.energy_level_beginning = pm.Param(
-                model.flexible_charging_points_set, initialize=0, mutable=True
-            )
-        else:
-            model.energy_level_beginning = pm.Param(
-                model.flexible_charging_points_set,
-                initialize=energy_level_beginning,
-                mutable=True,
-            )
-        model.FinalEVEnergyLevelEnd = pm.Constraint(
-            model.flexible_charging_points_set, model.time_end, rule=final_energy_level
-        )
-        model.FinalEVChargingPower = pm.Constraint(
-            model.flexible_charging_points_set,
-            model.time_end,
-            rule=final_charging_power,
-        )
-        if energy_level_end is None:
-            model.FinalEVEnergyLevelFix.deactivate()
-            model.FinalEVEnergyLevelEnd.deactivate()
-            model.FinalEVChargingPower.deactivate()
-        else:
-            if type(energy_level_end) != bool:
-                model.FinalEVEnergyLevelFix.deactivate()
-            elif type(energy_level_end) == bool:
-                model.FinalEVEnergyLevelEnd.deactivate()
-        # set initial charging power
-        charging_initial = kwargs.get("charging_start", None)
-        model.charging_initial = pm.Param(
-            model.flexible_charging_points_set,
-            initialize=charging_initial,
-            mutable=True,
-        )
-        model.slack_initial_charging_pos = pm.Var(
-            model.flexible_charging_points_set, bounds=(0, None)
-        )
-        model.slack_initial_charging_neg = pm.Var(
-            model.flexible_charging_points_set, bounds=(0, None)
-        )
-        model.InitialEVChargingPower = pm.Constraint(
-            model.flexible_charging_points_set,
-            model.time_zero,
-            rule=initial_charging_power,
-        )
-        if charging_initial is None:
-            model.InitialEVChargingPower.deactivate()
 
     if objective == "minimize_energy_level" or objective == "maximize_energy_level":
         model.AggrGrid = pm.Constraint(model.time_set, rule=aggregated_power)
@@ -884,20 +384,858 @@ def init_power_factors(model, branch, time):
     if kwargs.get("print_model", False):
         model.pprint()
     print("Successfully set up optimisation model.")
-    print("It took {} seconds to set up model.".format(perf_counter() - t1))
+    print(f"It took {perf_counter() - t1} seconds to set up model.")
     return model
 
 
-def set_lower_band_ev(model, cp, time):
-    return model.lower_ev_energy.loc[model.timeindex[time], cp]
-
-
-def set_upper_band_ev(model, cp, time):
-    return model.upper_ev_energy.loc[model.timeindex[time], cp]
+def add_grid_model_lopf(
+    model,
+    timeinvariant_parameters,
+    timesteps,
+    edisgo_object,
+    grid_object,
+    slack,
+    v_min,
+    v_max,
+    thermal_limits,
+    v_slack,
+    load_factor_rings,
+):
+    """
+    Method to add sets variables and constraints for including a representation of the
+    grid with a linearised power flow under omission of losses. Only applicable to
+    radial networks.
+    # Todo: add docstrings
+    """
 
+    def init_active_nodal_power(model, bus, time):
+        return (
+            timeinvariant_parameters["nodal_active_power"]
+            .T.loc[model.timeindex[time]]
+            .loc[bus]
+        )
 
-def set_power_band_ev(model, cp, time):
-    return model.upper_ev_power.loc[model.timeindex[time], cp]
+    def init_reactive_nodal_power(model, bus, time):
+        return (
+            timeinvariant_parameters["nodal_reactive_power"]
+            .T.loc[model.timeindex[time]]
+            .loc[bus]
+        )
+
+    def init_active_nodal_load(model, bus, time):
+        return (
+            timeinvariant_parameters["nodal_active_load"]
+            .T.loc[model.timeindex[time]]
+            .loc[bus]
+        )
+
+    def init_reactive_nodal_load(model, bus, time):
+        return (
+            timeinvariant_parameters["nodal_reactive_load"]
+            .T.loc[model.timeindex[time]]
+            .loc[bus]
+        )
+
+    def init_active_nodal_feedin(model, bus, time):
+        return (
+            timeinvariant_parameters["nodal_active_feedin"]
+            .T.loc[model.timeindex[time]]
+            .loc[bus]
+        )
+
+    def init_reactive_nodal_feedin(model, bus, time):
+        return (
+            timeinvariant_parameters["nodal_reactive_feedin"]
+            .T.loc[model.timeindex[time]]
+            .loc[bus]
+        )
+
+    def init_power_factors(model, branch, time):
+        return timeinvariant_parameters["power_factors"].loc[
+            branch, model.timeindex[time]
+        ]
+
+    # check if multiple voltage levels are present
+    if len(grid_object.buses_df.v_nom.unique()) > 1:
+        print(
+            "More than one voltage level included. Please make sure to "
+            "adapt all impedance values to one reference system."
+        )
+    # Sets and parameters
+    model.bus_set = pm.Set(initialize=grid_object.buses_df.index)
+    model.slack_bus = pm.Set(initialize=slack)
+    model.v_min = v_min
+    model.v_max = v_max
+    model.v_nom = timeinvariant_parameters["v_nom"]
+    model.thermal_limit = thermal_limits
+    model.pars = timeinvariant_parameters["pars"]
+    model.grid = grid_object
+    model.downstream_nodes_matrix = timeinvariant_parameters["downstream_nodes_matrix"]
+    model.nodal_active_power = pm.Param(
+        model.bus_set, model.time_set, initialize=init_active_nodal_power, mutable=True
+    )
+    model.nodal_reactive_power = pm.Param(
+        model.bus_set,
+        model.time_set,
+        initialize=init_reactive_nodal_power,
+        mutable=True,
+    )
+    model.nodal_active_load = pm.Param(
+        model.bus_set, model.time_set, initialize=init_active_nodal_load, mutable=True
+    )
+    model.nodal_reactive_load = pm.Param(
+        model.bus_set, model.time_set, initialize=init_reactive_nodal_load, mutable=True
+    )
+    model.nodal_active_feedin = pm.Param(
+        model.bus_set, model.time_set, initialize=init_active_nodal_feedin, mutable=True
+    )
+    model.nodal_reactive_feedin = pm.Param(
+        model.bus_set,
+        model.time_set,
+        initialize=init_reactive_nodal_feedin,
+        mutable=True,
+    )
+    model.tan_phi_load = timeinvariant_parameters["tan_phi_load"]
+    model.tan_phi_feedin = timeinvariant_parameters["tan_phi_feedin"]
+    model.v_slack = v_slack
+    model.branches = timeinvariant_parameters["branches"]
+    model.branch_set = pm.Set(initialize=model.branches.index)
+    model.underlying_branch_elements = timeinvariant_parameters[
+        "underlying_branch_elements"
+    ]
+    model.power_factors = pm.Param(
+        model.branch_set, model.time_set, initialize=init_power_factors, mutable=True
+    )
+    # add n-1 security # Todo: make optional?
+    # adapt i_lines_allowed for radial feeders
+    buses_in_cycles = list(
+        set(itertools.chain.from_iterable(edisgo_object.topology.rings))
+    )
+
+    # Find lines in cycles
+    lines_in_cycles = list(
+        grid_object.lines_df.loc[
+            grid_object.lines_df[["bus0", "bus1"]].isin(buses_in_cycles).all(axis=1)
+        ].index.values
+    )
+
+    model.branches_load_factors = pd.DataFrame(
+        index=model.time_set, columns=model.branch_set
+    )
+    model.branches_load_factors.loc[:, :] = 1
+    tmp_residual_load = edisgo_object.timeseries.residual_load.loc[timesteps]
+    indices = pd.DataFrame(index=timesteps, columns=["index"])
+    indices["index"] = [i for i in range(len(timesteps))]
+    model.branches_load_factors.loc[
+        indices.loc[tmp_residual_load.loc[timesteps] < 0].values.T[0], lines_in_cycles
+    ] = load_factor_rings  # Todo: distinction of mv and lv?
+    # Note: So far LV does not contain rings
+    # Variables
+    model.p_cum = pm.Var(model.branch_set, model.time_set)
+    model.slack_p_cum_pos = pm.Var(model.branch_set, model.time_set, bounds=(0, None))
+    model.slack_p_cum_neg = pm.Var(model.branch_set, model.time_set, bounds=(0, None))
+    model.q_cum = pm.Var(model.branch_set, model.time_set)
+    model.v = pm.Var(model.bus_set, model.time_set)
+    model.slack_v_pos = pm.Var(model.bus_set, model.time_set, bounds=(0, None))
+    model.slack_v_neg = pm.Var(model.bus_set, model.time_set, bounds=(0, None))
+    model.curtailment_load = pm.Var(
+        model.bus_set,
+        model.time_set,
+        bounds=lambda m, b, t: (0, m.nodal_active_load[b, t]),
+    )
+    model.curtailment_feedin = pm.Var(
+        model.bus_set,
+        model.time_set,
+        bounds=lambda m, b, t: (0, m.nodal_active_feedin[b, t]),
+    )
+    # add curtailment of flexible units if present
+    if hasattr(model, "flexible_charging_points_set"):
+        model.curtailment_ev = pm.Var(model.bus_set, model.time_set, bounds=(0, None))
+
+        model.UpperCurtEV = pm.Constraint(
+            model.bus_set, model.time_set, rule=upper_bound_curtailment_ev
+        )
+    # Constraints
+    model.ActivePower = pm.Constraint(
+        model.branch_set, model.time_set, rule=active_power
+    )
+    model.UpperActive = pm.Constraint(
+        model.branch_set, model.time_set, rule=upper_active_power
+    )
+    model.LowerActive = pm.Constraint(
+        model.branch_set, model.time_set, rule=lower_active_power
+    )
+    # model.ReactivePower = pm.Constraint(model.branch_set, model.time_set,
+    #                                     rule=reactive_power)
+    model.SlackVoltage = pm.Constraint(
+        model.slack_bus, model.time_set, rule=slack_voltage
+    )
+    model.VoltageDrop = pm.Constraint(
+        model.branch_set, model.time_set, rule=voltage_drop
+    )
+    model.UpperVoltage = pm.Constraint(
+        model.bus_set, model.time_set, rule=upper_voltage
+    )
+    model.LowerVoltage = pm.Constraint(
+        model.bus_set, model.time_set, rule=lower_voltage
+    )
+    # model.UpperCurtLoad = pm.Constraint(model.bus_set, model.time_set,
+    #                                     rule=upper_bound_curtailment_load)
+    return model
+
+
+def add_ev_model_bands(
+    model,
+    timeinvariant_parameters,
+    grid_object,
+    charging_efficiency,
+    energy_level_start=None,
+    energy_level_end=None,
+    energy_level_beginning=None,
+    charging_start=None,
+):
+    """
+    Method to add sets, variables and constraints for including EV flexibility in terms
+    of energy bands.
+    Todo: add docstrings
+    """
+    # Sets and parameters
+    model.charging_points_set = pm.Set(initialize=grid_object.charging_points_df.index)
+    model.flexible_charging_points_set = pm.Set(
+        initialize=timeinvariant_parameters["optimized_charging_points"]
+    )
+    model.inflexible_charging_points_set = (
+        model.charging_points_set - model.flexible_charging_points_set
+    )
+    model.upper_ev_power = timeinvariant_parameters["ev_flex_bands"]["upper_power"]
+    model.upper_ev_energy = timeinvariant_parameters["ev_flex_bands"]["upper_energy"]
+    model.lower_ev_energy = timeinvariant_parameters["ev_flex_bands"]["lower_energy"]
+    model.charging_efficiency = charging_efficiency
+    model.lower_bound_ev = pm.Param(
+        model.flexible_charging_points_set,
+        model.time_set,
+        initialize=set_lower_band_ev,
+        mutable=True,
+    )
+    model.upper_bound_ev = pm.Param(
+        model.flexible_charging_points_set,
+        model.time_set,
+        initialize=set_upper_band_ev,
+        mutable=True,
+    )
+    model.power_bound_ev = pm.Param(
+        model.flexible_charging_points_set,
+        model.time_set,
+        initialize=set_power_band_ev,
+        mutable=True,
+    )
+    # Variables
+    model.charging_ev = pm.Var(
+        model.flexible_charging_points_set,
+        model.time_set,
+        bounds=lambda m, b, t: (0, m.power_bound_ev[b, t]),
+    )
+
+    model.energy_level_ev = pm.Var(
+        model.flexible_charging_points_set,
+        model.time_set,
+        bounds=lambda m, b, t: (m.lower_bound_ev[b, t], m.upper_bound_ev[b, t]),
+    )
+    # Constraints
+    model.EVCharging = pm.Constraint(
+        model.flexible_charging_points_set, model.time_non_zero, rule=charging_ev
+    )
+    # set initial energy level, Todo: possibly add own method for rolling horizon
+    model.energy_level_start = pm.Param(
+        model.flexible_charging_points_set,
+        initialize=energy_level_start,
+        mutable=True,
+        within=pm.Any,
+    )
+    model.slack_initial_energy_pos = pm.Var(
+        model.flexible_charging_points_set, bounds=(0, None)
+    )
+    model.slack_initial_energy_neg = pm.Var(
+        model.flexible_charging_points_set, bounds=(0, None)
+    )
+    model.InitialEVEnergyLevel = pm.Constraint(
+        model.flexible_charging_points_set,
+        model.time_zero,
+        rule=initial_energy_level,
+    )
+    model.InitialEVEnergyLevelStart = pm.Constraint(
+        model.flexible_charging_points_set, model.time_zero, rule=fixed_energy_level
+    )
+    if energy_level_start is None:
+        model.InitialEVEnergyLevel.deactivate()
+    else:
+        model.InitialEVEnergyLevelStart.deactivate()
+    # set final energy level and if necessary charging power
+    model.energy_level_end = pm.Param(
+        model.flexible_charging_points_set,
+        initialize=energy_level_end,
+        mutable=True,
+        within=pm.Any,
+    )
+    model.FinalEVEnergyLevelFix = pm.Constraint(
+        model.flexible_charging_points_set, model.time_end, rule=fixed_energy_level
+    )
+
+    if energy_level_beginning is None:
+        model.energy_level_beginning = pm.Param(
+            model.flexible_charging_points_set, initialize=0, mutable=True
+        )
+    else:
+        model.energy_level_beginning = pm.Param(
+            model.flexible_charging_points_set,
+            initialize=energy_level_beginning,
+            mutable=True,
+        )
+    model.FinalEVEnergyLevelEnd = pm.Constraint(
+        model.flexible_charging_points_set, model.time_end, rule=final_energy_level
+    )
+    model.FinalEVChargingPower = pm.Constraint(
+        model.flexible_charging_points_set,
+        model.time_end,
+        rule=final_charging_power,
+    )
+    if energy_level_end is None:
+        model.FinalEVEnergyLevelFix.deactivate()
+        model.FinalEVEnergyLevelEnd.deactivate()
+        model.FinalEVChargingPower.deactivate()
+    else:
+        if type(energy_level_end) != bool:
+            model.FinalEVEnergyLevelFix.deactivate()
+        elif type(energy_level_end) == bool:
+            model.FinalEVEnergyLevelEnd.deactivate()
+    # set initial charging power
+    model.charging_initial = pm.Param(
+        model.flexible_charging_points_set,
+        initialize=charging_start,
+        mutable=True,
+    )
+    model.slack_initial_charging_pos = pm.Var(
+        model.flexible_charging_points_set, bounds=(0, None)
+    )
+    model.slack_initial_charging_neg = pm.Var(
+        model.flexible_charging_points_set, bounds=(0, None)
+    )
+    model.InitialEVChargingPower = pm.Constraint(
+        model.flexible_charging_points_set,
+        model.time_zero,
+        rule=initial_charging_power,
+    )
+    if charging_start is None:
+        model.InitialEVChargingPower.deactivate()
+
+    return model
+
+
+def update_model(
+    model, timesteps, parameters, optimize_storage=True, optimize_ev=True, **kwargs
+):
+    """
+    Method to update model parameter where necessary if rolling horizon
+    optimization is chosen.
+
+    Parameters
+    ----------
+    model
+    timesteps
+    parameters
+    optimize_storage
+    optimize_ev
+    kwargs
+
+    Returns
+    -------
+
+    """
+    print("Updating model")
+    t1 = perf_counter()
+    for i in model.time_set:
+        overlap = i - len(timesteps) + 1
+        if overlap > 0:
+            timeindex = timesteps[-1] + pd.to_timedelta(model.time_increment) * overlap
+            indexer = timesteps[-1]
+        else:
+            timeindex = timesteps[i]
+            indexer = timesteps[i]
+        model.timeindex[i].set_value(timeindex)
+        model.residual_load[i].set_value(
+            parameters["res_load_inflexible_units"][indexer]
+        )
+        for bus in model.bus_set:
+            model.nodal_active_power[bus, i].set_value(
+                parameters["nodal_active_power"].loc[bus, indexer]
+            )
+            model.nodal_reactive_power[bus, i].set_value(
+                parameters["nodal_reactive_power"].loc[bus, indexer]
+            )
+            model.nodal_active_load[bus, i].set_value(
+                parameters["nodal_active_load"].loc[bus, indexer]
+            )
+            model.nodal_reactive_load[bus, i].set_value(
+                parameters["nodal_reactive_load"].loc[bus, indexer]
+            )
+            model.nodal_active_feedin[bus, i].set_value(
+                parameters["nodal_active_feedin"].loc[bus, indexer]
+            )
+            model.nodal_reactive_feedin[bus, i].set_value(
+                parameters["nodal_reactive_feedin"].loc[bus, indexer]
+            )
+
+        for branch in model.branch_set:
+            model.power_factors[branch, i].set_value(
+                parameters["power_factors"].loc[branch, indexer]
+            )
+
+    if optimize_ev:
+        for t in model.time_set:
+            overlap = t - len(timesteps) + 1
+            if overlap > 0:
+                indexer = len(timesteps) - 1
+            else:
+                indexer = t
+            for cp in model.flexible_charging_points_set:
+                model.power_bound_ev[cp, t].set_value(
+                    set_power_band_ev(model, cp, indexer)
+                )
+                model.lower_bound_ev[cp, t].set_value(
+                    set_lower_band_ev(model, cp, indexer)
+                )
+                model.upper_bound_ev[cp, t].set_value(
+                    set_upper_band_ev(model, cp, indexer)
+                )
+        # set initial energy level
+        energy_level_start = kwargs.get("energy_level_start", None)
+        charging_initial = kwargs.get("charging_start", None)
+        # if run is new start of era deactivate initial energy level,
+        # otherwise activate initial energy and charging
+        if energy_level_start is None:
+            model.InitialEVEnergyLevel.deactivate()
+            model.InitialEVEnergyLevelStart.activate()
+        else:
+            for cp in model.flexible_charging_points_set:
+                model.energy_level_start[cp].set_value(energy_level_start[cp])
+            model.InitialEVEnergyLevel.activate()
+            model.InitialEVEnergyLevelStart.deactivate()
+        # set initial charging
+        if charging_initial is not None:
+            for cp in model.flexible_charging_points_set:
+                model.charging_initial[cp].set_value(charging_initial[cp])
+            model.InitialEVChargingPower.activate()
+        # set energy level beginning if necessary
+
+        energy_level_beginning = kwargs.get("energy_level_beginning", None)
+        if energy_level_beginning is not None:
+            for cp in model.flexible_charging_points_set:
+                model.energy_level_beginning[cp].set_value(energy_level_beginning[cp])
+
+        # set final energy level and if necessary charging power
+        energy_level_end = kwargs.get("energy_level_end", None)
+        if energy_level_end is None:
+            model.FinalEVEnergyLevelFix.deactivate()
+            model.FinalEVEnergyLevelEnd.deactivate()
+            model.FinalEVChargingPower.deactivate()
+        elif type(energy_level_end) == bool:
+            model.FinalEVEnergyLevelFix.activate()
+            model.FinalEVEnergyLevelEnd.deactivate()
+            model.FinalEVChargingPower.activate()
+        else:
+            for cp in model.flexible_charging_points_set:
+                model.energy_level_end[cp].set_value(energy_level_end[cp])
+            model.FinalEVEnergyLevelEnd.activate()
+            model.FinalEVChargingPower.activate()
+            model.FinalEVEnergyLevelFix.deactivate()
+
+    if optimize_storage:
+        raise NotImplementedError
+    print("It took {} seconds to update the model.".format(perf_counter() - t1))
+    return model
+
+
+def optimize(model, solver, load_solutions=True, mode=None):
+    """
+    Method to run the optimization and extract the results.
+
+    :param model: pyomo.environ.ConcreteModel
+    :param solver: str
+                    Solver type, e.g. 'glpk', 'gurobi', 'ipopt'
+    :param save_dir: str
+                    directory to which results are saved, default None will
+                    no saving of the results
+    :return:
+    """
+    print("Starting optimisation")
+    t1 = perf_counter()
+    opt = pm.SolverFactory(solver)
+    opt.options["threads"] = 16
+
+    # Optimize
+    results = opt.solve(model, tee=True, load_solutions=load_solutions)
+
+    if (results.solver.status == SolverStatus.ok) and (
+        results.solver.termination_condition == TerminationCondition.optimal
+    ):
+        print("Model Solved to Optimality")
+        # Extract results
+        time_dict = {t: model.timeindex[t].value for t in model.time_set}
+        result_dict = {}
+        if hasattr(model, "storage_set"):
+            result_dict["x_charge"] = (
+                pd.Series(model.charging.extract_values())
+                .unstack()
+                .rename(columns=time_dict)
+                .T
+            )
+            result_dict["soc"] = (
+                pd.Series(model.soc.extract_values())
+                .unstack()
+                .rename(columns=time_dict)
+                .T
+            )
+        if hasattr(model, "flexible_charging_points_set"):
+            result_dict["x_charge_ev"] = (
+                pd.Series(model.charging_ev.extract_values())
+                .unstack()
+                .rename(columns=time_dict)
+                .T
+            )
+            result_dict["energy_level_cp"] = (
+                pd.Series(model.energy_level_ev.extract_values())
+                .unstack()
+                .rename(columns=time_dict)
+                .T
+            )
+            result_dict["slack_charging"] = pd.Series(
+                model.slack_initial_charging_pos.extract_values()
+            ) + pd.Series(model.slack_initial_charging_neg.extract_values())
+            result_dict["slack_energy"] = pd.Series(
+                model.slack_initial_energy_pos.extract_values()
+            ) + pd.Series(model.slack_initial_energy_neg.extract_values())
+        result_dict["curtailment_load"] = (
+            pd.Series(model.curtailment_load.extract_values())
+            .unstack()
+            .rename(columns=time_dict)
+            .T
+        )
+        result_dict["curtailment_feedin"] = (
+            pd.Series(model.curtailment_feedin.extract_values())
+            .unstack()
+            .rename(columns=time_dict)
+            .T
+        )
+        result_dict["curtailment_ev"] = (
+            pd.Series(model.curtailment_ev.extract_values())
+            .unstack()
+            .rename(columns=time_dict)
+            .T
+        )
+        result_dict["p_line"] = (
+            pd.Series(model.p_cum.extract_values())
+            .unstack()
+            .rename(columns=time_dict)
+            .T
+        )
+        result_dict["q_line"] = (
+            pd.Series(model.q_cum.extract_values())
+            .unstack()
+            .rename(columns=time_dict)
+            .T
+        )
+        result_dict["v_bus"] = (
+            pd.Series(model.v.extract_values())
+            .unstack()
+            .rename(columns=time_dict)
+            .T.apply(np.sqrt)
+        )
+        result_dict["slack_v_pos"] = (
+            pd.Series(model.slack_v_pos.extract_values())
+            .unstack()
+            .rename(columns=time_dict)
+            .T
+        )
+        result_dict["slack_v_neg"] = (
+            pd.Series(model.slack_v_neg.extract_values())
+            .unstack()
+            .rename(columns=time_dict)
+            .T
+        )
+        result_dict["slack_p_cum_pos"] = (
+            pd.Series(model.slack_p_cum_pos.extract_values())
+            .unstack()
+            .rename(columns=time_dict)
+            .T
+        )
+        result_dict["slack_p_cum_neg"] = (
+            pd.Series(model.slack_p_cum_pos.extract_values())
+            .unstack()
+            .rename(columns=time_dict)
+            .T
+        )
+        if mode == "energy_band":
+            result_dict["p_aggr"] = pd.Series(
+                model.grid_power_flexible.extract_values()
+            ).rename(time_dict)
+        # Todo: check if this works
+        index = result_dict["curtailment_load"].index[
+            result_dict["curtailment_load"].index.isin(model.tan_phi_load.index)
+        ]
+        result_dict["curtailment_reactive_load"] = (
+            result_dict["curtailment_load"]
+            .multiply(
+                model.tan_phi_load.loc[index, result_dict["curtailment_load"].columns]
+            )
+            .dropna(how="all")
+        )
+        result_dict["curtailment_reactive_feedin"] = (
+            result_dict["curtailment_feedin"]
+            .multiply(
+                model.tan_phi_feedin.loc[
+                    index, result_dict["curtailment_feedin"].columns
+                ]
+            )
+            .dropna(how="all")
+        )
+
+        print("It took {} seconds to optimize model.".format(perf_counter() - t1))
+        return result_dict
+    elif results.solver.termination_condition == TerminationCondition.infeasible:
+        print("Model is infeasible")
+        return
+        # Do something when model in infeasible
+    else:
+        print("Solver Status: ", results.solver.status)
+        return
+
+
+def setup_grid_object(object):
+    """
+    Set up the grid and edisgo object.
+    """
+    if hasattr(object, "topology"):
+        grid_object = deepcopy(object.topology)
+        edisgo_object = deepcopy(object)
+        slack = grid_object.mv_grid.station.index
+    else:
+        grid_object = deepcopy(object)
+        edisgo_object = deepcopy(object.edisgo_obj)
+        slack = [
+            grid_object.transformers_df.bus1.iloc[0]
+        ]  # Todo: careful with MV grid, does not work with that right?
+    return edisgo_object, grid_object, slack
+
+
+def concat_parallel_branch_elements(grid_object):
+    """
+    Method to merge parallel lines and transformers into one element, respectively.
+
+    Parameters
+    ----------
+    grid_object
+
+    Returns
+    -------
+
+    """
+    lines = fuse_parallel_branches(grid_object.lines_df)
+    trafos = grid_object.transformers_df.loc[
+        grid_object.transformers_df.bus0.isin(grid_object.buses_df.index)
+    ].loc[grid_object.transformers_df.bus1.isin(grid_object.buses_df.index)]
+    transformers = fuse_parallel_branches(trafos)
+    return pd.concat([lines, transformers], sort=False)
+
+
+def fuse_parallel_branches(branches):
+    branches_tmp = branches[["bus0", "bus1"]]
+    parallel_branches = pd.DataFrame(columns=branches.columns)
+    if branches_tmp.duplicated().any():
+        duplicated_branches = branches_tmp.loc[branches_tmp.duplicated(keep=False)]
+        duplicated_branches["visited"] = False
+        branches_tmp.drop(duplicated_branches.index, inplace=True)
+        for name, buses in duplicated_branches.iterrows():
+            if duplicated_branches.loc[name, "visited"]:
+                continue
+            else:
+                parallel_branches_tmp = duplicated_branches.loc[
+                    (duplicated_branches == buses).all(axis=1)
+                ]
+                duplicated_branches.loc[parallel_branches_tmp.index, "visited"] = True
+                name_par = "_".join(str.split(name, "_")[:-1])
+                parallel_branches.loc[name_par] = branches.loc[name]
+                parallel_branches.loc[
+                    name_par, ["r", "x", "s_nom"]
+                ] = calculate_impedance_for_parallel_components(
+                    branches.loc[parallel_branches_tmp.index, ["r", "x", "s_nom"]],
+                    pu=False,
+                )
+    fused_branches = pd.concat(
+        [branches.loc[branches_tmp.index], parallel_branches], sort=False
+    )
+    return fused_branches
+
+
+def get_underlying_elements(parameters):
+    def _get_underlying_elements(
+        downstream_elements, power_factors, parameters, branch
+    ):
+        bus0 = parameters["branches"].loc[branch, "bus0"]
+        bus1 = parameters["branches"].loc[branch, "bus1"]
+        s_nom = parameters["branches"].loc[branch, parameters["pars"]["s_nom"]]
+        relevant_buses_bus0 = (
+            parameters["downstream_nodes_matrix"]
+            .loc[bus0][parameters["downstream_nodes_matrix"].loc[bus0] == 1]
+            .index.values
+        )
+        relevant_buses_bus1 = (
+            parameters["downstream_nodes_matrix"]
+            .loc[bus1][parameters["downstream_nodes_matrix"].loc[bus1] == 1]
+            .index.values
+        )
+        relevant_buses = list(
+            set(relevant_buses_bus0).intersection(relevant_buses_bus1)
+        )
+        downstream_elements.loc[branch, "buses"] = relevant_buses
+        if (
+            parameters["nodal_reactive_power"]
+            .loc[relevant_buses]
+            .sum()
+            .divide(s_nom)
+            .apply(abs)
+            > 1
+        ).any():
+            print(
+                "Careful: Reactive power already exceeding line capacity for branch "
+                "{}.".format(branch)
+            )
+        power_factors.loc[branch] = (
+            1
+            - parameters["nodal_reactive_power"]
+            .loc[relevant_buses]
+            .sum()
+            .divide(s_nom)
+            .apply(np.square)
+        ).apply(np.sqrt)
+        if "optimized_storage_units" in parameters:
+            downstream_elements.loc[branch, "flexible_storage"] = (
+                parameters["grid_object"]
+                .storage_units_df.loc[
+                    parameters["grid_object"].storage_units_df.index.isin(
+                        parameters["optimized_storage_units"]
+                    )
+                    & parameters["grid_object"].storage_units_df.bus.isin(
+                        relevant_buses
+                    )
+                ]
+                .index.values
+            )
+        else:
+            downstream_elements.loc[branch, "flexible_storage"] = []
+        if "optimized_charging_points" in parameters:
+            downstream_elements.loc[branch, "flexible_ev"] = (
+                parameters["grid_object"]
+                .charging_points_df.loc[
+                    parameters["grid_object"].charging_points_df.index.isin(
+                        parameters["optimized_charging_points"]
+                    )
+                    & parameters["grid_object"].charging_points_df.bus.isin(
+                        relevant_buses
+                    )
+                ]
+                .index.values
+            )
+        else:
+            downstream_elements.loc[branch, "flexible_ev"] = []
+        return downstream_elements, power_factors
+
+    downstream_elements = pd.DataFrame(
+        index=parameters["branches"].index,
+        columns=["buses", "flexible_storage", "flexible_ev"],
+    )
+    power_factors = pd.DataFrame(
+        index=parameters["branches"].index,
+        columns=parameters["nodal_active_power"].columns,
+    )
+    for branch in downstream_elements.index:
+        downstream_elements, power_factors = _get_underlying_elements(
+            downstream_elements, power_factors, parameters, branch
+        )
+    if power_factors.isna().any().any():
+        print(
+            "WARNING: Branch {} is overloaded with reactive power. Still needs "
+            "handling.".format(branch)
+        )
+        power_factors = power_factors.fillna(
+            0.01
+        )  # Todo: ask Gaby and Birgit about this
+    return downstream_elements, power_factors
+
+
+def get_residual_load_of_not_optimized_components(
+    grid,
+    edisgo,
+    relevant_storage_units=None,
+    relevant_charging_points=None,
+    relevant_generators=None,
+    relevant_loads=None,
+):
+    """
+    Method to get residual load of fixed components.
+
+    Parameters
+    ----------
+    grid
+    edisgo
+    relevant_storage_units
+    relevant_charging_points
+    relevant_generators
+    relevant_loads
+
+    Returns
+    -------
+
+    """
+    if relevant_loads is None:
+        relevant_loads = grid.loads_df.index
+    if relevant_generators is None:
+        relevant_generators = grid.generators_df.index
+    if relevant_storage_units is None:
+        relevant_storage_units = grid.storage_units_df.index
+    if relevant_charging_points is None:
+        relevant_charging_points = grid.charging_points_df.index
+
+    if edisgo.timeseries.charging_points_active_power.empty:
+        return (
+            edisgo.timeseries.generators_active_power[relevant_generators].sum(axis=1)
+            + edisgo.timeseries.storage_units_active_power[relevant_storage_units].sum(
+                axis=1
+            )
+            - edisgo.timeseries.loads_active_power[relevant_loads].sum(axis=1)
+        ).loc[edisgo.timeseries.timeindex]
+    else:
+        return (
+            edisgo.timeseries.generators_active_power[relevant_generators].sum(axis=1)
+            + edisgo.timeseries.storage_units_active_power[relevant_storage_units].sum(
+                axis=1
+            )
+            - edisgo.timeseries.loads_active_power[relevant_loads].sum(axis=1)
+            - edisgo.timeseries.charging_points_active_power[
+                relevant_charging_points
+            ].sum(axis=1)
+        ).loc[edisgo.timeseries.timeindex]
+
+
+def set_lower_band_ev(model, cp, time):
+    return model.lower_ev_energy.loc[model.timeindex[time], cp]
+
+
+def set_upper_band_ev(model, cp, time):
+    return model.upper_ev_energy.loc[model.timeindex[time], cp]
+
+
+def set_power_band_ev(model, cp, time):
+    return model.upper_ev_power.loc[model.timeindex[time], cp]
 
 
 def active_power(model, branch, time):
@@ -1418,410 +1756,131 @@ def minimize_residual_load(model):
                 1e-2
                 * (
                     model.curtailment_load[bus, time]
-                    + model.curtailment_feedin[bus, time]
-                    + 0.5 * model.curtailment_ev[bus, time]
-                )
-                + 1000 * (model.slack_v_pos[bus, time] + model.slack_v_neg[bus, time])
-                for bus in model.bus_set
-                for time in model.time_set
-            )
-            + 1000 * (slack_charging + slack_energy)
-            + 1000
-            * sum(
-                model.slack_p_cum_pos[branch, time]
-                + model.slack_p_cum_neg[branch, time]
-                for branch in model.branch_set
-                for time in model.time_set
-            )
-        )
-    else:
-        return (
-            1e-5 * sum(model.grid_residual_load[time] ** 2 for time in model.time_set)
-            + sum(
-                1e-2
-                * (
-                    model.curtailment_load[bus, time]
-                    + model.curtailment_feedin[bus, time]
-                )
-                + 1000 * (model.slack_v_pos[bus, time] + model.slack_v_neg[bus, time])
-                for bus in model.bus_set
-                for time in model.time_set
-            )
-            + 1000
-            * sum(
-                model.slack_p_cum_pos[branch, time]
-                + model.slack_p_cum_neg[branch, time]
-                for branch in model.branch_set
-                for time in model.time_set
-            )
-        )
-
-
-def minimize_loading(model):
-    """
-    Objective minimizing curtailment and squared term of component loading
-    :param model:
-    :return:
-    """
-    if hasattr(model, "slack_initial_charging_pos"):
-        slack_charging = sum(
-            model.slack_initial_charging_pos[cp] + model.slack_initial_charging_neg[cp]
-            for cp in model.flexible_charging_points_set
-        )
-    else:
-        slack_charging = 0
-    if hasattr(model, "slack_initial_energy_pos"):
-        slack_energy = sum(
-            model.slack_initial_energy_pos[cp] + model.slack_initial_energy_neg[cp]
-            for cp in model.flexible_charging_points_set
-        )
-    else:
-        slack_energy = 0
-    if hasattr(model, "charging_points_set"):
-        return (
-            1e-5
-            * sum(
-                (
-                    model.p_cum[branch, time]
-                    / (
-                        model.power_factors[branch, time]
-                        * model.branches.loc[branch, model.pars["s_nom"]]
-                    )
-                )
-                ** 2
-                for branch in model.branch_set
-                for time in model.time_set
-            )
-            + sum(
-                1e-2
-                * (
-                    model.curtailment_load[bus, time]
-                    + model.curtailment_feedin[bus, time]
-                    + 0.5 * model.curtailment_ev[bus, time]
-                )
-                + 1000 * (model.slack_v_pos[bus, time] + model.slack_v_neg[bus, time])
-                for bus in model.bus_set
-                for time in model.time_set
-            )
-            + 1000 * (slack_charging + slack_energy)
-            + 1000
-            * sum(
-                model.slack_p_cum_pos[branch, time]
-                + model.slack_p_cum_neg[branch, time]
-                for branch in model.branch_set
-                for time in model.time_set
-            )
-        )
-    else:
-        return (
-            1e-5
-            * sum(
-                model.p_cum[branch, time].divide(
-                    model.power_factors[branch, time]
-                    * model.branches.loc[branch, model.pars["s_nom"]]
-                )
-                ** 2
-                for branch in model.branch_set
-                for time in model.time_set
-            )
-            + sum(
-                1e-2
-                * (
-                    model.curtailment_load[bus, time]
-                    + model.curtailment_feedin[bus, time]
-                )
-                + 1000 * (model.slack_v_pos[bus, time] + model.slack_v_neg[bus, time])
-                for bus in model.bus_set
-                for time in model.time_set
-            )
-            + 1000
-            * sum(
-                model.slack_p_cum_pos[branch, time]
-                + model.slack_p_cum_neg[branch, time]
-                for branch in model.branch_set
-                for time in model.time_set
-            )
-        )
-
-
-def update_model(
-    model, timesteps, parameters, optimize_storage=True, optimize_ev=True, **kwargs
-):
-    """
-    Method to update model parameter where necessary if rolling horizon
-    optimization is chosen.
-
-    Parameters
-    ----------
-    model
-    timesteps
-    parameters
-    optimize_storage
-    optimize_ev
-    kwargs
-
-    Returns
-    -------
-
-    """
-    print("Updating model")
-    t1 = perf_counter()
-    for i in model.time_set:
-        overlap = i - len(timesteps) + 1
-        if overlap > 0:
-            timeindex = timesteps[-1] + pd.to_timedelta(model.time_increment) * overlap
-            indexer = timesteps[-1]
-        else:
-            timeindex = timesteps[i]
-            indexer = timesteps[i]
-        model.timeindex[i].set_value(timeindex)
-        model.residual_load[i].set_value(
-            parameters["res_load_inflexible_units"][indexer]
-        )
-        for bus in model.bus_set:
-            model.nodal_active_power[bus, i].set_value(
-                parameters["nodal_active_power"].loc[bus, indexer]
-            )
-            model.nodal_reactive_power[bus, i].set_value(
-                parameters["nodal_reactive_power"].loc[bus, indexer]
-            )
-            model.nodal_active_load[bus, i].set_value(
-                parameters["nodal_active_load"].loc[bus, indexer]
-            )
-            model.nodal_reactive_load[bus, i].set_value(
-                parameters["nodal_reactive_load"].loc[bus, indexer]
-            )
-            model.nodal_active_feedin[bus, i].set_value(
-                parameters["nodal_active_feedin"].loc[bus, indexer]
-            )
-            model.nodal_reactive_feedin[bus, i].set_value(
-                parameters["nodal_reactive_feedin"].loc[bus, indexer]
-            )
-
-        for branch in model.branch_set:
-            model.power_factors[branch, i].set_value(
-                parameters["power_factors"].loc[branch, indexer]
-            )
-
-    if optimize_ev:
-        for t in model.time_set:
-            overlap = t - len(timesteps) + 1
-            if overlap > 0:
-                indexer = len(timesteps) - 1
-            else:
-                indexer = t
-            for cp in model.flexible_charging_points_set:
-                model.power_bound_ev[cp, t].set_value(
-                    set_power_band_ev(model, cp, indexer)
-                )
-                model.lower_bound_ev[cp, t].set_value(
-                    set_lower_band_ev(model, cp, indexer)
-                )
-                model.upper_bound_ev[cp, t].set_value(
-                    set_upper_band_ev(model, cp, indexer)
-                )
-        # set initial energy level
-        energy_level_start = kwargs.get("energy_level_start", None)
-        charging_initial = kwargs.get("charging_start", None)
-        # if run is new start of era deactivate initial energy level,
-        # otherwise activate initial energy and charging
-        if energy_level_start is None:
-            model.InitialEVEnergyLevel.deactivate()
-            model.InitialEVEnergyLevelStart.activate()
-        else:
-            for cp in model.flexible_charging_points_set:
-                model.energy_level_start[cp].set_value(energy_level_start[cp])
-            model.InitialEVEnergyLevel.activate()
-            model.InitialEVEnergyLevelStart.deactivate()
-        # set initial charging
-        if charging_initial is not None:
-            for cp in model.flexible_charging_points_set:
-                model.charging_initial[cp].set_value(charging_initial[cp])
-            model.InitialEVChargingPower.activate()
-        # set energy level beginning if necessary
-
-        energy_level_beginning = kwargs.get("energy_level_beginning", None)
-        if energy_level_beginning is not None:
-            for cp in model.flexible_charging_points_set:
-                model.energy_level_beginning[cp].set_value(energy_level_beginning[cp])
-
-        # set final energy level and if necessary charging power
-        energy_level_end = kwargs.get("energy_level_end", None)
-        if energy_level_end is None:
-            model.FinalEVEnergyLevelFix.deactivate()
-            model.FinalEVEnergyLevelEnd.deactivate()
-            model.FinalEVChargingPower.deactivate()
-        elif type(energy_level_end) == bool:
-            model.FinalEVEnergyLevelFix.activate()
-            model.FinalEVEnergyLevelEnd.deactivate()
-            model.FinalEVChargingPower.activate()
-        else:
-            for cp in model.flexible_charging_points_set:
-                model.energy_level_end[cp].set_value(energy_level_end[cp])
-            model.FinalEVEnergyLevelEnd.activate()
-            model.FinalEVChargingPower.activate()
-            model.FinalEVEnergyLevelFix.deactivate()
-
-    if optimize_storage:
-        raise NotImplementedError
-    print("It took {} seconds to update the model.".format(perf_counter() - t1))
-    return model
+                    + model.curtailment_feedin[bus, time]
+                    + 0.5 * model.curtailment_ev[bus, time]
+                )
+                + 1000 * (model.slack_v_pos[bus, time] + model.slack_v_neg[bus, time])
+                for bus in model.bus_set
+                for time in model.time_set
+            )
+            + 1000 * (slack_charging + slack_energy)
+            + 1000
+            * sum(
+                model.slack_p_cum_pos[branch, time]
+                + model.slack_p_cum_neg[branch, time]
+                for branch in model.branch_set
+                for time in model.time_set
+            )
+        )
+    else:
+        return (
+            1e-5 * sum(model.grid_residual_load[time] ** 2 for time in model.time_set)
+            + sum(
+                1e-2
+                * (
+                    model.curtailment_load[bus, time]
+                    + model.curtailment_feedin[bus, time]
+                )
+                + 1000 * (model.slack_v_pos[bus, time] + model.slack_v_neg[bus, time])
+                for bus in model.bus_set
+                for time in model.time_set
+            )
+            + 1000
+            * sum(
+                model.slack_p_cum_pos[branch, time]
+                + model.slack_p_cum_neg[branch, time]
+                for branch in model.branch_set
+                for time in model.time_set
+            )
+        )
 
 
-def optimize(model, solver, load_solutions=True, mode=None):
+def minimize_loading(model):
     """
-    Method to run the optimization and extract the results.
-
-    :param model: pyomo.environ.ConcreteModel
-    :param solver: str
-                    Solver type, e.g. 'glpk', 'gurobi', 'ipopt'
-    :param save_dir: str
-                    directory to which results are saved, default None will
-                    no saving of the results
+    Objective minimizing curtailment and squared term of component loading
+    :param model:
     :return:
     """
-    print("Starting optimisation")
-    t1 = perf_counter()
-    opt = pm.SolverFactory(solver)
-    opt.options["threads"] = 16
-
-    # Optimize
-    results = opt.solve(model, tee=True, load_solutions=load_solutions)
-
-    if (results.solver.status == SolverStatus.ok) and (
-        results.solver.termination_condition == TerminationCondition.optimal
-    ):
-        print("Model Solved to Optimality")
-        # Extract results
-        time_dict = {t: model.timeindex[t].value for t in model.time_set}
-        result_dict = {}
-        if hasattr(model, "storage_set"):
-            result_dict["x_charge"] = (
-                pd.Series(model.charging.extract_values())
-                .unstack()
-                .rename(columns=time_dict)
-                .T
-            )
-            result_dict["soc"] = (
-                pd.Series(model.soc.extract_values())
-                .unstack()
-                .rename(columns=time_dict)
-                .T
+    if hasattr(model, "slack_initial_charging_pos"):
+        slack_charging = sum(
+            model.slack_initial_charging_pos[cp] + model.slack_initial_charging_neg[cp]
+            for cp in model.flexible_charging_points_set
+        )
+    else:
+        slack_charging = 0
+    if hasattr(model, "slack_initial_energy_pos"):
+        slack_energy = sum(
+            model.slack_initial_energy_pos[cp] + model.slack_initial_energy_neg[cp]
+            for cp in model.flexible_charging_points_set
+        )
+    else:
+        slack_energy = 0
+    if hasattr(model, "charging_points_set"):
+        return (
+            1e-5
+            * sum(
+                (
+                    model.p_cum[branch, time]
+                    / (
+                        model.power_factors[branch, time]
+                        * model.branches.loc[branch, model.pars["s_nom"]]
+                    )
+                )
+                ** 2
+                for branch in model.branch_set
+                for time in model.time_set
             )
-        if hasattr(model, "flexible_charging_points_set"):
-            result_dict["x_charge_ev"] = (
-                pd.Series(model.charging_ev.extract_values())
-                .unstack()
-                .rename(columns=time_dict)
-                .T
+            + sum(
+                1e-2
+                * (
+                    model.curtailment_load[bus, time]
+                    + model.curtailment_feedin[bus, time]
+                    + 0.5 * model.curtailment_ev[bus, time]
+                )
+                + 1000 * (model.slack_v_pos[bus, time] + model.slack_v_neg[bus, time])
+                for bus in model.bus_set
+                for time in model.time_set
             )
-            result_dict["energy_level_cp"] = (
-                pd.Series(model.energy_level_ev.extract_values())
-                .unstack()
-                .rename(columns=time_dict)
-                .T
+            + 1000 * (slack_charging + slack_energy)
+            + 1000
+            * sum(
+                model.slack_p_cum_pos[branch, time]
+                + model.slack_p_cum_neg[branch, time]
+                for branch in model.branch_set
+                for time in model.time_set
             )
-            result_dict["slack_charging"] = pd.Series(
-                model.slack_initial_charging_pos.extract_values()
-            ) + pd.Series(model.slack_initial_charging_neg.extract_values())
-            result_dict["slack_energy"] = pd.Series(
-                model.slack_initial_energy_pos.extract_values()
-            ) + pd.Series(model.slack_initial_energy_neg.extract_values())
-        result_dict["curtailment_load"] = (
-            pd.Series(model.curtailment_load.extract_values())
-            .unstack()
-            .rename(columns=time_dict)
-            .T
-        )
-        result_dict["curtailment_feedin"] = (
-            pd.Series(model.curtailment_feedin.extract_values())
-            .unstack()
-            .rename(columns=time_dict)
-            .T
-        )
-        result_dict["curtailment_ev"] = (
-            pd.Series(model.curtailment_ev.extract_values())
-            .unstack()
-            .rename(columns=time_dict)
-            .T
-        )
-        result_dict["p_line"] = (
-            pd.Series(model.p_cum.extract_values())
-            .unstack()
-            .rename(columns=time_dict)
-            .T
-        )
-        result_dict["q_line"] = (
-            pd.Series(model.q_cum.extract_values())
-            .unstack()
-            .rename(columns=time_dict)
-            .T
-        )
-        result_dict["v_bus"] = (
-            pd.Series(model.v.extract_values())
-            .unstack()
-            .rename(columns=time_dict)
-            .T.apply(np.sqrt)
-        )
-        result_dict["slack_v_pos"] = (
-            pd.Series(model.slack_v_pos.extract_values())
-            .unstack()
-            .rename(columns=time_dict)
-            .T
-        )
-        result_dict["slack_v_neg"] = (
-            pd.Series(model.slack_v_neg.extract_values())
-            .unstack()
-            .rename(columns=time_dict)
-            .T
-        )
-        result_dict["slack_p_cum_pos"] = (
-            pd.Series(model.slack_p_cum_pos.extract_values())
-            .unstack()
-            .rename(columns=time_dict)
-            .T
-        )
-        result_dict["slack_p_cum_neg"] = (
-            pd.Series(model.slack_p_cum_pos.extract_values())
-            .unstack()
-            .rename(columns=time_dict)
-            .T
         )
-        if mode == "energy_band":
-            result_dict["p_aggr"] = pd.Series(
-                model.grid_power_flexible.extract_values()
-            ).rename(time_dict)
-        # Todo: check if this works
-        index = result_dict["curtailment_load"].index[
-            result_dict["curtailment_load"].index.isin(model.tan_phi_load.index)
-        ]
-        result_dict["curtailment_reactive_load"] = (
-            result_dict["curtailment_load"]
-            .multiply(
-                model.tan_phi_load.loc[index, result_dict["curtailment_load"].columns]
+    else:
+        return (
+            1e-5
+            * sum(
+                model.p_cum[branch, time].divide(
+                    model.power_factors[branch, time]
+                    * model.branches.loc[branch, model.pars["s_nom"]]
+                )
+                ** 2
+                for branch in model.branch_set
+                for time in model.time_set
             )
-            .dropna(how="all")
-        )
-        result_dict["curtailment_reactive_feedin"] = (
-            result_dict["curtailment_feedin"]
-            .multiply(
-                model.tan_phi_feedin.loc[
-                    index, result_dict["curtailment_feedin"].columns
-                ]
+            + sum(
+                1e-2
+                * (
+                    model.curtailment_load[bus, time]
+                    + model.curtailment_feedin[bus, time]
+                )
+                + 1000 * (model.slack_v_pos[bus, time] + model.slack_v_neg[bus, time])
+                for bus in model.bus_set
+                for time in model.time_set
+            )
+            + 1000
+            * sum(
+                model.slack_p_cum_pos[branch, time]
+                + model.slack_p_cum_neg[branch, time]
+                for branch in model.branch_set
+                for time in model.time_set
             )
-            .dropna(how="all")
         )
 
-        print("It took {} seconds to optimize model.".format(perf_counter() - t1))
-        return result_dict
-    elif results.solver.termination_condition == TerminationCondition.infeasible:
-        print("Model is infeasible")
-        return
-        # Do something when model in infeasible
-    else:
-        print("Solver Status: ", results.solver.status)
-        return
-
 
 def combine_results_for_grid(feeders, grid_id, res_dir, res_name):
     res_grid = pd.DataFrame()
@@ -1850,5 +1909,3 @@ def combine_results_for_grid(feeders, grid_id, res_dir, res_name):
             print("Feeder {} not added".format(feeder_id))
     res_grid = res_grid.loc[~res_grid.index.duplicated(keep="last")]
     return res_grid
-
-# fmt: on
diff --git a/edisgo/opf/timeseries_reduction.py b/edisgo/opf/timeseries_reduction.py
index ede525bf2..31013f244 100644
--- a/edisgo/opf/timeseries_reduction.py
+++ b/edisgo/opf/timeseries_reduction.py
@@ -1,20 +1,19 @@
 import logging
-import pandas as pd
+
 import numpy as np
+import pandas as pd
+from sklearn.cluster import KMeans
+from sklearn.preprocessing import StandardScaler
 
 from edisgo.flex_opt import check_tech_constraints
-
-from sklearn.preprocessing import StandardScaler
-from sklearn.cluster import KMeans
+from results_helper_functions import relative_load, voltage_diff
 
 logger = logging.getLogger(__name__)
 
 
 def _scored_critical_current(edisgo_obj, grid):
     # Get allowed current per line per time step
-    i_lines_allowed = check_tech_constraints.lines_allowed_load(
-        edisgo_obj, "mv"
-    )
+    i_lines_allowed = check_tech_constraints.lines_allowed_load(edisgo_obj, "mv")
     i_lines_pfa = edisgo_obj.results.i_res[grid.lines_df.index]
 
     # Get current relative to allowed current
@@ -32,6 +31,29 @@ def _scored_critical_current(edisgo_obj, grid):
     return crit_lines_score.sort_values(ascending=False)
 
 
+def _scored_most_critical_current(edisgo_obj):
+
+    # Get current relative to allowed current
+    relative_i_res = relative_load(edisgo_obj)
+
+    # Get only timesteps that are maximum loading of at least one component
+    # relative_i_res_tmp = relative_i_res.loc[list(set(relative_i_res.idxmax().values))]
+
+    # Get lines that have violations
+    crit_lines_score = relative_i_res[relative_i_res > 1]
+
+    # Remove time steps with no violations
+    crit_lines_score = crit_lines_score.dropna(how="all", axis=0)
+
+    # Get most critical timesteps per component
+    if not crit_lines_score.empty:
+        idx_max = list(set(crit_lines_score.idxmax().dropna().values))
+
+        crit_lines_score = crit_lines_score.loc[idx_max] - 1
+    crit_lines_score = crit_lines_score.sum(axis=1)
+    return crit_lines_score.sort_values(ascending=False)
+
+
 def _scored_critical_overvoltage(edisgo_obj, grid):
     nodes = grid.buses_df.index
 
@@ -43,38 +65,88 @@ def _scored_critical_overvoltage(edisgo_obj, grid):
         edisgo_obj, voltage_levels="mv"
     )
     _, voltage_diff_ov = check_tech_constraints.voltage_diff(
-        edisgo_obj,
-        nodes,
-        v_dev_allowed_upper,
-        v_dev_allowed_lower
+        edisgo_obj, nodes, v_dev_allowed_upper, v_dev_allowed_lower
     )
 
     # Get score for nodes that are over or under the allowed deviations
     voltage_diff_ov = (
-        voltage_diff_ov[voltage_diff_ov > 0]
-        .dropna(axis=1, how="all")
-        .sum(axis=0)
+        voltage_diff_ov[voltage_diff_ov > 0].dropna(axis=1, how="all").sum(axis=0)
     )
     return voltage_diff_ov.sort_values(ascending=False)
 
 
-def get_steps_curtailment(edisgo_obj, percentage=0.5):
+def _scored_most_critical_voltage_issues(edisgo_obj):
+    voltage_diff_ov = voltage_diff(edisgo_obj)
+
+    # Get score for nodes that are over or under the allowed deviations
+    voltage_diff_ov = voltage_diff_ov.abs()[voltage_diff_ov.abs() > 0].dropna(
+        axis=0, how="all"
+    )
+    if not voltage_diff_ov.empty:
+        # Get index of most critical timesteps
+        idx_max = list(set(voltage_diff_ov.idxmax().dropna().values))
+        # Get sum of all violations at most critical steps
+        voltage_diff_ov = voltage_diff_ov.loc[idx_max]
+    voltage_diff_ov = voltage_diff_ov.sum(axis=1)
+
+    return voltage_diff_ov.sort_values(ascending=False)
+
+
+def get_steps_reinforcement(edisgo_obj, percentage=1.0):
     """
     Get the time steps with the most critical violations for curtailment
     optimization.
-
     Parameters
     -----------
     edisgo_obj : :class:`~.EDisGo`
         The eDisGo API object
     percentage : float
         The percentage of most critical time steps to select
-
     Returns
     --------
     `pandas.DatetimeIndex`
         the reduced time index for modeling curtailment
+    """
+    # Run power flow if not available
+    if edisgo_obj.results.i_res is None or edisgo_obj.results.i_res.empty:
+        logger.debug("Running initial power flow")
+        edisgo_obj.analyze()
 
+    # Select most critical steps based on current violations
+    current_scores = _scored_most_critical_current(edisgo_obj)
+    num_steps_current = int(len(current_scores) * percentage)
+    steps = current_scores[:num_steps_current].index.tolist()
+
+    # Select most critical steps based on voltage violations
+    voltage_scores = _scored_most_critical_voltage_issues(edisgo_obj)
+    num_steps_voltage = int(len(voltage_scores) * percentage)
+    steps.extend(
+        voltage_scores[:num_steps_voltage].index.tolist()
+    )  # Todo: Can this cause duplicated?
+
+    if len(steps) == 0:
+        logger.warning("No critical steps detected. No network expansion required.")
+
+    # Strip duplicates
+    steps = list(dict.fromkeys(steps))
+
+    return pd.DatetimeIndex(steps)
+
+
+def get_steps_curtailment(edisgo_obj, percentage=0.5):
+    """
+    Get the time steps with the most critical violations for curtailment
+    optimization.
+    Parameters
+    -----------
+    edisgo_obj : :class:`~.EDisGo`
+        The eDisGo API object
+    percentage : float
+        The percentage of most critical time steps to select
+    Returns
+    --------
+    `pandas.DatetimeIndex`
+        the reduced time index for modeling curtailment
     """
 
     # Run power flow if not available
@@ -98,9 +170,7 @@ def get_steps_curtailment(edisgo_obj, percentage=0.5):
     steps.extend(get_steps_storage(edisgo_obj, window=0).tolist())
 
     if len(steps) == 0:
-        logger.warning(
-            "No critical steps detected. No network expansion required."
-        )
+        logger.warning("No critical steps detected. No network expansion required.")
 
     # Strip duplicates
     steps = list(dict.fromkeys(steps))
@@ -111,7 +181,6 @@ def get_steps_curtailment(edisgo_obj, percentage=0.5):
 def get_steps_storage(edisgo_obj, window=5):
     """
     Get the most critical time steps from series for storage problems.
-
     Parameters
     -----------
     edisgo_obj : :class:`~.EDisGo`
@@ -119,12 +188,10 @@ def get_steps_storage(edisgo_obj, window=5):
     window : int
         The additional hours to include before and after each critical time
         step.
-
     Returns
     -------
     `pandas.DatetimeIndex`
         the reduced time index for modeling storage
-
     """
     # Run power flow if not available
     if edisgo_obj.results.i_res is None:
@@ -141,14 +208,14 @@ def get_steps_storage(edisgo_obj, window=5):
         nodes = pd.DataFrame(v)
         if "time_index" in nodes:
             for step in nodes["time_index"]:
-                if not step in crit_periods:
+                if step not in crit_periods:
                     crit_periods.append(step)
 
     # Get periods with current violations
     crit_lines = check_tech_constraints.mv_line_load(edisgo_obj)
     if "time_index" in crit_lines:
         for step in crit_lines["time_index"]:
-            if not step in crit_periods:
+            if step not in crit_periods:
                 crit_periods.append(step)
 
     reduced = []
@@ -164,9 +231,7 @@ def get_steps_storage(edisgo_obj, window=5):
     reduced = list(dict.fromkeys(reduced))
 
     if len(reduced) == 0:
-        logger.warning(
-            "No critical steps detected. No network expansion required."
-        )
+        logger.warning("No critical steps detected. No network expansion required.")
 
     return pd.DatetimeIndex(reduced)
 
@@ -175,7 +240,6 @@ def get_linked_steps(cluster_params, num_steps=24, keep_steps=[]):
     """
     Use provided data to identify representative time steps and create mapping
     Dict that can be passed to optimization
-
     Parameters
     -----------
     cluster_params : :pandas:`pandas.DataFrame<DataFrame>`
@@ -186,13 +250,11 @@ def get_linked_steps(cluster_params, num_steps=24, keep_steps=[]):
     keep_steps : Iterable of the same type as cluster_params.index
         Time steps to retain with full resolution, regardless of
         clustering result.
-
     Returns
     -------
     dict
         Dictionary where each represented time step is a key and its
         representative time step is a value.
-
     """
 
     # From all values, find the subvector with the smallest SSD to a given
diff --git a/edisgo/tools/plots.py b/edisgo/tools/plots.py
index 92ffe0609..287b2a75e 100644
--- a/edisgo/tools/plots.py
+++ b/edisgo/tools/plots.py
@@ -155,7 +155,7 @@ def add_basemap(ax, zoom=12):
     Adds map to a plot.
 
     """
-    url = ctx.sources.ST_TONER_LITE
+    url = ctx.providers.Stamen.TonerLite
     xmin, xmax, ymin, ymax = ax.axis()
     basemap, extent = ctx.bounds2img(
         xmin, ymin, xmax, ymax, zoom=zoom, source=url
diff --git a/examples/example_grid_reinforcement.py b/examples/example_grid_reinforcement.py
index b7c326fea..fa7d8f425 100644
--- a/examples/example_grid_reinforcement.py
+++ b/examples/example_grid_reinforcement.py
@@ -25,33 +25,41 @@
 
 """
 
+import logging
 import os
+
 import pandas as pd
 import requests
 
 from edisgo import EDisGo
 from edisgo.network.results import Results
 
-import logging
-logger = logging.getLogger('edisgo')
+logger = logging.getLogger("edisgo")
 logger.setLevel(logging.DEBUG)
 
 
 def run_example():
     # Specify path to directory containing ding0 grid csv files
-    edisgo_path = os.path.join(os.path.expanduser('~'), '.eDisGo')
-    dingo_grid_path = os.path.join(edisgo_path,
-                                   'ding0_example_grid')
+    edisgo_path = os.path.join(os.path.expanduser("~"), ".eDisGo")
+    dingo_grid_path = os.path.join(edisgo_path, "ding0_example_grid")
     # Download example grid data in case it does not yet exist
     if not os.path.isdir(dingo_grid_path):
         logger.debug("Download example grid data.")
         os.makedirs(dingo_grid_path)
-        file_list = ["buses.csv", "lines.csv",
-                     "generators.csv", "loads.csv",
-                     "transformers.csv", "transformers_hvmv.csv",
-                     "switches.csv", "network.csv"]
-        base_path = "https://raw.githubusercontent.com/openego/eDisGo/" \
-                    "dev/tests/ding0_test_network_2/"
+        file_list = [
+            "buses.csv",
+            "lines.csv",
+            "generators.csv",
+            "loads.csv",
+            "transformers.csv",
+            "transformers_hvmv.csv",
+            "switches.csv",
+            "network.csv",
+        ]
+        base_path = (
+            "https://raw.githubusercontent.com/openego/eDisGo/"
+            "dev/tests/ding0_test_network_2/"
+        )
         for f in file_list:
             file = os.path.join(dingo_grid_path, f)
             req = requests.get(base_path + f)
@@ -59,29 +67,26 @@ def run_example():
                 fout.write(req.content)
 
     # Set scenario to define future power plant capacities
-    scenario = 'nep2035'
+    scenario = "nep2035"
 
     # Set up worst-case scenario
-    edisgo = EDisGo(ding0_grid=dingo_grid_path,
-                    worst_case_analysis='worst-case')
+    edisgo = EDisGo(ding0_grid=dingo_grid_path, worst_case_analysis="worst-case")
 
     # Reinforce ding0 grid to obtain a stable status quo grid
-    logging.info('Conduct grid reinforcement to obtain stable '
-                 'status quo grid.')
+    logging.info("Conduct grid reinforcement to obtain stable " "status quo grid.")
     # Overwrite config parameters for allowed voltage deviations in
     # initial topology reinforcement to better represent currently used limits
-    edisgo.config[
-        'grid_expansion_allowed_voltage_deviations'] = {
-        'feedin_case_lower': 0.9,
-        'load_case_upper': 1.1,
-        'hv_mv_trafo_offset': 0.04,
-        'hv_mv_trafo_control_deviation': 0.0,
-        'mv_load_case_max_v_deviation': 0.055,
-        'mv_feedin_case_max_v_deviation': 0.02,
-        'lv_load_case_max_v_deviation': 0.065,
-        'lv_feedin_case_max_v_deviation': 0.03,
-        'mv_lv_station_load_case_max_v_deviation': 0.02,
-        'mv_lv_station_feedin_case_max_v_deviation': 0.01
+    edisgo.config["grid_expansion_allowed_voltage_deviations"] = {
+        "feed-in_case_lower": 0.9,
+        "load_case_upper": 1.1,
+        "hv_mv_trafo_offset": 0.04,
+        "hv_mv_trafo_control_deviation": 0.0,
+        "mv_load_case_max_v_deviation": 0.055,
+        "mv_feed-in_case_max_v_deviation": 0.02,
+        "lv_load_case_max_v_deviation": 0.065,
+        "lv_feed-in_case_max_v_deviation": 0.03,
+        "mv_lv_station_load_case_max_v_deviation": 0.02,
+        "mv_lv_station_feed-in_case_max_v_deviation": 0.01,
     }
     # Conduct reinforcement
     edisgo.reinforce()
@@ -91,7 +96,7 @@ def run_example():
     edisgo.config = None
 
     # Calculate expansion costs for NEP scenario
-    logging.info('Determine grid expansion costs for NEP scenario.')
+    logging.info("Determine grid expansion costs for NEP scenario.")
 
     # Get data on generators in NEP scenario and connect generators to the grid
     edisgo.import_generators(generator_scenario=scenario)
@@ -102,20 +107,18 @@ def run_example():
     # Get total grid expansion costs
     total_costs = edisgo.results.grid_expansion_costs.total_costs.sum()
     logging.info(
-        'Grid expansion costs for NEP scenario are: {} kEUR.'.format(
-            total_costs)
+        "Grid expansion costs for NEP scenario are: {} kEUR.".format(total_costs)
     )
 
     # Save grid expansion results
     edisgo.results.to_csv(
-        directory=dingo_grid_path,
-        parameters={'grid_expansion_results': None}
+        directory=dingo_grid_path, parameters={"grid_expansion_results": None}
     )
 
-    logging.info('SUCCESS!')
+    logging.info("SUCCESS!")
 
     return total_costs
 
 
-if __name__ == '__main__':
+if __name__ == "__main__":
     run_example()
diff --git a/examples/example_linopt.py b/examples/example_linopt.py
new file mode 100644
index 000000000..da4cf881f
--- /dev/null
+++ b/examples/example_linopt.py
@@ -0,0 +1,77 @@
+import os
+
+import matplotlib.pyplot as plt
+import pandas as pd
+
+import edisgo.opf.lopf as opt
+from edisgo.edisgo import import_edisgo_from_files
+from edisgo.network.electromobility import (Electromobility,
+                                            get_energy_bands_for_optimization)
+from edisgo.network.timeseries import get_component_timeseries
+from Script_prepare_grids_for_optimization import \
+    get_downstream_nodes_matrix_iterative
+
+grid_dir = "minimum_working"
+opt_ev = True
+opt_stor = False
+save_res = False
+
+if os.path.isfile("x_charge_ev_pre.csv"):
+    ts_pre = pd.read_csv("x_charge_ev_pre.csv", index_col=0, parse_dates=True)
+else:
+    ts_pre = pd.DataFrame()
+
+timeindex = pd.date_range("2011-01-01", periods=8760, freq="h")
+storage_ts = pd.DataFrame({"Storage 1": 8760 * [0]}, index=timeindex)
+
+edisgo = import_edisgo_from_files(grid_dir)
+get_component_timeseries(
+    edisgo,
+    timeseries_load="demandlib",
+    timeseries_generation_fluctuating="oedb",
+    timeseries_storage_units=storage_ts,
+)
+timesteps = edisgo.timeseries.timeindex[7 * 24 : 2 * 24 * 7]
+
+cp_id = 1
+ev_data = pd.read_csv(
+    os.path.join(grid_dir, "BEV_standing_times_minimum_working.csv"), index_col=0
+)
+charging_events = ev_data.loc[ev_data.chargingdemand > 0]
+charging_events["charging_park_id"] = cp_id
+Electromobility(edisgo_obj=edisgo)
+edisgo.electromobility.charging_processes_df = charging_events
+cp = edisgo.add_component(
+    "ChargingPoint", bus="Bus 2", p_nom=0.011, use_case="home", add_ts=False
+)
+edisgo.electromobility.integrated_charging_parks_df = pd.DataFrame(
+    index=[cp_id], columns=["edisgo_id"], data=cp
+)
+edisgo.electromobility.simbev_config_df = pd.DataFrame(
+    index=["eta_CP", "stepsize"], columns=["value"], data=[0.9, 60]
+)
+energy_bands = get_energy_bands_for_optimization(edisgo_obj=edisgo, use_case="home")
+
+downstream_node_matrix = get_downstream_nodes_matrix_iterative(edisgo.topology)
+parameters = opt.prepare_time_invariant_parameters(
+    edisgo,
+    downstream_node_matrix,
+    pu=False,
+    optimize_storage=False,
+    optimize_ev_charging=True,
+    ev_flex_bands=energy_bands,
+)
+model = opt.setup_model(
+    parameters, timesteps=timesteps, objective="residual_load", optimize_storage=False
+)
+
+results = opt.optimize(model, "gurobi")
+results["x_charge_ev"].plot()
+plt.show()
+if not ts_pre.empty:
+    ts_pre.plot()
+    plt.show()
+    pd.testing.assert_frame_equal(ts_pre, results["x_charge_ev"])
+if save_res:
+    results["x_charge_ev"].to_csv("x_charge_ev_pre.csv")
+print("SUCCESS")
diff --git a/results_helper_functions.py b/results_helper_functions.py
index c2adbf9f6..ef7887de2 100644
--- a/results_helper_functions.py
+++ b/results_helper_functions.py
@@ -1,7 +1,10 @@
-import pandas as pd
 import numpy as np
-from edisgo.flex_opt.check_tech_constraints import lines_allowed_load, \
-    lines_relative_load, _mv_allowed_voltage_limits, _lv_allowed_voltage_limits
+import pandas as pd
+
+from edisgo.flex_opt.check_tech_constraints import (_lv_allowed_voltage_limits,
+                                                    _mv_allowed_voltage_limits,
+                                                    lines_allowed_load,
+                                                    lines_relative_load)
 
 
 def relative_load(edisgo_obj):
@@ -36,11 +39,19 @@ def relative_load(edisgo_obj):
     if any(lv_lines.isin(edisgo_obj.results.i_res.columns)):
         lv_lines_allowed_load = lines_allowed_load(edisgo_obj, "lv")
         allowed_load_lines = pd.concat(
-            [allowed_load_lines,
-             lv_lines_allowed_load.loc[:, edisgo_obj.results.i_res.columns[
-                                              edisgo_obj.results.i_res.columns.isin(
-                                                  lv_lines_allowed_load.columns)]]],
-            axis=1)
+            [
+                allowed_load_lines,
+                lv_lines_allowed_load.loc[
+                    :,
+                    edisgo_obj.results.i_res.columns[
+                        edisgo_obj.results.i_res.columns.isin(
+                            lv_lines_allowed_load.columns
+                        )
+                    ],
+                ],
+            ],
+            axis=1,
+        )
 
     # calculated relative load for lines
     rel_load = lines_relative_load(edisgo_obj, allowed_load_lines)
@@ -49,7 +60,8 @@ def relative_load(edisgo_obj):
     # check if power flow was conducted for stations
     if not edisgo_obj.results.s_res.empty:
         if not edisgo_obj.results.s_res.loc[
-           :, edisgo_obj.results.s_res.columns.str.contains("Transformer")].empty:
+            :, edisgo_obj.results.s_res.columns.str.contains("Transformer")
+        ].empty:
 
             load_factor = edisgo_obj.timeseries.timesteps_load_feedin_case.apply(
                 lambda _: edisgo_obj.config["grid_expansion_load_factors"][
@@ -69,8 +81,9 @@ def relative_load(edisgo_obj):
                     # get maximum allowed apparent power of station in each time
                     # step
                     s_station_allowed = sum(transformers_df.s_nom) * load_factor
-                    rel_load["mvlv_station_{}".format(grid)] = \
+                    rel_load["mvlv_station_{}".format(grid)] = (
                         s_station_pfa / s_station_allowed
+                    )
 
     # HV-MV station
     # check if power flow was conducted for MV
@@ -88,8 +101,9 @@ def relative_load(edisgo_obj):
             ]
         )
         s_station_allowed = s_station * load_factor
-        rel_load["hvmv_station_{}".format(edisgo_obj.topology.mv_grid)] = \
+        rel_load["hvmv_station_{}".format(edisgo_obj.topology.mv_grid)] = (
             s_station_pfa / s_station_allowed
+        )
 
     return rel_load
 
@@ -106,23 +120,23 @@ def voltage_diff_stations(edisgo_obj):
     voltage_base = edisgo_obj.results.v_res.loc[:, primary_sides.values]
 
     v_allowed_per_case = {}
-    v_allowed_per_case["feedin_case_upper"] = (
-            voltage_base
-            + edisgo_obj.config["grid_expansion_allowed_voltage_deviations"][
-                "mv_lv_station_feedin_case_max_v_deviation"
-            ]
+    v_allowed_per_case["feed-in_case_upper"] = (
+        voltage_base
+        + edisgo_obj.config["grid_expansion_allowed_voltage_deviations"][
+            "mv_lv_station_feed-in_case_max_v_deviation"
+        ]
     )
     v_allowed_per_case["load_case_lower"] = (
-            voltage_base
-            - edisgo_obj.config["grid_expansion_allowed_voltage_deviations"][
-                "mv_lv_station_load_case_max_v_deviation"
-            ]
+        voltage_base
+        - edisgo_obj.config["grid_expansion_allowed_voltage_deviations"][
+            "mv_lv_station_load_case_max_v_deviation"
+        ]
     )
 
     timeindex = voltage_base.index
-    v_allowed_per_case["feedin_case_lower"] = pd.DataFrame(
+    v_allowed_per_case["feed-in_case_lower"] = pd.DataFrame(
         edisgo_obj.config["grid_expansion_allowed_voltage_deviations"][
-            "feedin_case_lower"
+            "feed-in_case_lower"
         ],
         columns=v_allowed_per_case["load_case_lower"].columns,
         index=timeindex,
@@ -148,28 +162,24 @@ def voltage_diff_stations(edisgo_obj):
         )
 
     # rename columns to secondary side
-    rename_dict = {primary_sides[g]: secondary_sides[g] for g in
-                   edisgo_obj.topology.mv_grid.lv_grids}
+    rename_dict = {
+        primary_sides[g]: secondary_sides[g]
+        for g in edisgo_obj.topology.mv_grid.lv_grids
+    }
     v_dev_allowed_upper_all.rename(columns=rename_dict, inplace=True)
     v_dev_allowed_lower_all.rename(columns=rename_dict, inplace=True)
 
-    v_mag_pu_pfa_all = edisgo_obj.results.v_res.loc[:,
-                       v_dev_allowed_upper_all.columns]
+    v_mag_pu_pfa_all = edisgo_obj.results.v_res.loc[:, v_dev_allowed_upper_all.columns]
 
-    overvoltage = v_mag_pu_pfa_all[
-        v_mag_pu_pfa_all > v_dev_allowed_upper_all
-        ]
-    undervoltage = v_mag_pu_pfa_all[
-        v_mag_pu_pfa_all < v_dev_allowed_lower_all
-        ]
+    overvoltage = v_mag_pu_pfa_all[v_mag_pu_pfa_all > v_dev_allowed_upper_all]
+    undervoltage = v_mag_pu_pfa_all[v_mag_pu_pfa_all < v_dev_allowed_lower_all]
 
     # overvoltage diff (positive)
     overvoltage_diff = overvoltage - v_dev_allowed_upper_all
     # undervoltage diff (negative)
     undervoltage_diff = undervoltage - v_dev_allowed_lower_all
 
-    voltage_difference_all = overvoltage_diff.fillna(
-        0) + undervoltage_diff.fillna(0)
+    voltage_difference_all = overvoltage_diff.fillna(0) + undervoltage_diff.fillna(0)
 
     return voltage_difference_all
 
@@ -201,10 +211,12 @@ def voltage_diff(edisgo_obj):
     mv_buses = edisgo_obj.topology.mv_grid.buses_df.index
     if any(mv_buses.isin(edisgo_obj.results.v_res.columns)):
         v_dev_allowed_upper, v_dev_allowed_lower = _mv_allowed_voltage_limits(
-            edisgo_obj, voltage_levels="mv")
+            edisgo_obj, voltage_levels="mv"
+        )
 
         v_mag_pu_pfa = edisgo_obj.results.v_res.loc[
-                       :, edisgo_obj.topology.mv_grid.buses_df.index]
+            :, edisgo_obj.topology.mv_grid.buses_df.index
+        ]
 
         v_dev_allowed_upper_format = np.tile(
             (v_dev_allowed_upper.loc[v_mag_pu_pfa.index]).values,
@@ -214,19 +226,14 @@ def voltage_diff(edisgo_obj):
             (v_dev_allowed_lower.loc[v_mag_pu_pfa.index]).values,
             (v_mag_pu_pfa.shape[1], 1),
         )
-        overvoltage = v_mag_pu_pfa.T[
-            v_mag_pu_pfa.T > v_dev_allowed_upper_format
-        ]
-        undervoltage = v_mag_pu_pfa.T[
-            v_mag_pu_pfa.T < v_dev_allowed_lower_format
-        ]
+        overvoltage = v_mag_pu_pfa.T[v_mag_pu_pfa.T > v_dev_allowed_upper_format]
+        undervoltage = v_mag_pu_pfa.T[v_mag_pu_pfa.T < v_dev_allowed_lower_format]
 
         # overvoltage diff (positive)
         overvoltage_diff = overvoltage - v_dev_allowed_upper_format
         # undervoltage diff (negative)
         undervoltage_diff = undervoltage - v_dev_allowed_lower_format
-        voltage_difference = (overvoltage_diff.fillna(0) +
-                              undervoltage_diff.fillna(0))
+        voltage_difference = overvoltage_diff.fillna(0) + undervoltage_diff.fillna(0)
         voltage_difference = voltage_difference.T
     else:
         voltage_difference = pd.DataFrame()
@@ -235,13 +242,12 @@ def voltage_diff(edisgo_obj):
     # check if power flow was conducted for stations
     if not edisgo_obj.results.s_res.empty:
         if not edisgo_obj.results.s_res.loc[
-               :,
-               edisgo_obj.results.s_res.columns.str.contains("Transformer")].empty:
+            :, edisgo_obj.results.s_res.columns.str.contains("Transformer")
+        ].empty:
 
             voltage_difference_stations = voltage_diff_stations(edisgo_obj)
             voltage_difference = pd.concat(
-                [voltage_difference, voltage_difference_stations],
-                sort=False, axis=1
+                [voltage_difference, voltage_difference_stations], sort=False, axis=1
             )
 
     # LV buses
@@ -254,15 +260,13 @@ def voltage_diff(edisgo_obj):
         for lv_grid in edisgo_obj.topology.mv_grid.lv_grids:
 
             # check if grid was included in power flow
-            if any(lv_grid.lines_df.index.isin(
-                    edisgo_obj.results.s_res.columns)):
+            if any(lv_grid.lines_df.index.isin(edisgo_obj.results.s_res.columns)):
 
-                v_dev_allowed_upper, v_dev_allowed_lower = \
-                    _lv_allowed_voltage_limits(
-                        edisgo_obj, lv_grid, mode=None)
+                v_dev_allowed_upper, v_dev_allowed_lower = _lv_allowed_voltage_limits(
+                    edisgo_obj, lv_grid, mode=None
+                )
 
-                v_mag_pu_pfa = edisgo_obj.results.v_res.loc[
-                               :, lv_grid.buses_df.index]
+                v_mag_pu_pfa = edisgo_obj.results.v_res.loc[:, lv_grid.buses_df.index]
 
                 v_dev_allowed_upper_format = np.tile(
                     (v_dev_allowed_upper.loc[v_mag_pu_pfa.index]).values,
@@ -274,23 +278,29 @@ def voltage_diff(edisgo_obj):
                 )
                 overvoltage = v_mag_pu_pfa.T[
                     v_mag_pu_pfa.T > v_dev_allowed_upper_format
-                    ]
+                ]
                 undervoltage = v_mag_pu_pfa.T[
                     v_mag_pu_pfa.T < v_dev_allowed_lower_format
-                    ]
+                ]
 
                 # overvoltage diff (positive)
                 overvoltage_diff = overvoltage - v_dev_allowed_upper_format
                 # undervoltage diff (negative)
                 undervoltage_diff = undervoltage - v_dev_allowed_lower_format
-                voltage_difference_lv_grid = (
-                        overvoltage_diff.fillna(0) +
-                        undervoltage_diff.fillna(0)
-                )
+                voltage_difference_lv_grid = overvoltage_diff.fillna(
+                    0
+                ) + undervoltage_diff.fillna(0)
                 voltage_difference = pd.concat(
-                    [voltage_difference,
-                     voltage_difference_lv_grid[~voltage_difference_lv_grid.T.columns.isin(voltage_difference.columns)].T],
-                    sort=False, axis=1
+                    [
+                        voltage_difference,
+                        voltage_difference_lv_grid[
+                            ~voltage_difference_lv_grid.T.columns.isin(
+                                voltage_difference.columns
+                            )
+                        ].T,
+                    ],
+                    sort=False,
+                    axis=1,
                 )
 
-    return voltage_difference
\ No newline at end of file
+    return voltage_difference