Seeking help with implementing cost function for flexible electricity price and solar storage in OEMOF model

Hello everyone,

I am currently working on an energy system model using OEMOF and facing the challenge of implementing a cost function for a flexible electricity price as well as integrating a solar storage system. Specifically, I aim to calculate the total cost of electricity over a certain period and examine the influence of a flexible electricity price from the supplier.

My objectives are:

  1. Implementing a cost function that accounts for the flexible electricity price.
  2. Comparing the total costs with and without solar storage over the same period.
  3. Analyzing the impact of a flexible electricity price from the supplier on the total costs.

I have already started developing the code but am stuck at some points. Here is a snippet of my current code:

###########################################################################
# imports
###########################################################################
import pandas as pd
import matplotlib.pyplot as plt
import base64
from io import BytesIO
import logging
import os
import pprint as pp
import warnings
from datetime import datetime
import pandas as pd

import matplotlib.pyplot as plt
import pandas as pd
from oemof.tools import logger
import oemof.solph as solph
import matplotlib.pyplot as plt
from oemof.solph import flows

from oemof.solph import EnergySystem
from oemof.solph import Model
from oemof.solph import buses
from oemof.solph import components as cmp
from oemof.solph import create_time_index
from oemof.solph import Flow
from oemof.solph import helpers
from oemof.solph import processing
from oemof.solph import views
import numpy as np
from oemof.solph import constraints

def main():
    # *************************************************************************
    # ********** PART 1 - Define and optimise the energy system ***************
    # *************************************************************************

    # Read data file

    verbrauchsdaten_mit_pv_file = os.path.join(os.getcwd(), 'verbrauchsdaten_mit_pv.csv')
    supplier_cost_file = os.path.join(os.getcwd(), 'Day_Night_price_Supplier_.csv')                                         

    verbrauchsdaten_mit_pv_data = pd.read_csv(verbrauchsdaten_mit_pv_file)
    supplier_cost_data = pd.read_csv(supplier_cost_file)
   
    solver = "cbc"  # 'glpk', 'gurobi',....
    debug = False  # Set number_of_timesteps to 3 to get a readable lp-file.
    solver_verbose = False  # show/hide solver output

    # initiate the logger (see the API docs for more information)
    logger.define_logging(
        logfile="oemof_example.log",
        screen_level=logging.INFO,
        file_level=logging.INFO,
    )

    date_time_index = pd.date_range('1/1/2024 00:00:00', periods=48, freq='30min')

    # create the energysystem and assign the time index
    energysystem = EnergySystem(
        timeindex=date_time_index, infer_last_interval=False 
    )

##########################################################################
    # Create oemof objects
    ##########################################################################

    logging.info("Create oemof objects")

    # The bus objects were assigned to variables which makes it easier to
    # connect components to these buses (see below).

    # create electricity bus
    bel = buses.Bus(label="electricity")

    # adding the buses to the energy system
    energysystem.add(bel)

    # create excess component for the electricity bus to allow overproduction
    energysystem.add(cmp.Sink(label="excess_bel", inputs={bel: flows.Flow()}))

    energysystem.add(
        cmp.Source(
            label="supplier",
            outputs={
                bel: flows.Flow(
                    nominal_value=1000, 
                    variable_costs=supplier_cost_data["costs"], 
                    custom_attributes={"emission_factor": 0.20}, 
                )
            },
        )
    )

    # create fixed source object representing pv power plants
    energysystem.add(
        cmp.Source(
            label="pv",
            outputs={bel: flows.Flow(fix=verbrauchsdaten_mit_pv_data["pv"], nominal_value=1000)},
        )                                                             
    )                                                                      

    # create simple sink object representing the electrical demand 1
    energysystem.add(
        cmp.Sink(
            label="consumer1",
            inputs={bel: flows.Flow(fix=verbrauchsdaten_mit_pv_data["consumer1"], nominal_value=1, custom_attributes={"emission_factor": 0.60},)},
        )                                                                 
    )

    # create simple sink object representing the electrical demand 2
    energysystem.add(
        cmp.Sink(
            label="consumer2",
            inputs={bel: flows.Flow(fix=verbrauchsdaten_mit_pv_data["consumer2"], nominal_value=1, custom_attributes={"emission_factor": 0.30},)} 
        )                                                                   
    )

    # create storage object representing a battery
    storage = cmp.GenericStorage(
        nominal_storage_capacity=25000,
        label="storage",
        inputs={bel: flows.Flow(nominal_value=25000 / 6)},
        outputs={
            bel: flows.Flow(nominal_value=25000 / 6, variable_costs=0.001)
        },
        loss_rate=0.00,
        initial_storage_level=None,
        inflow_conversion_factor=1,
        outflow_conversion_factor=0.8,
    )

    energysystem.add(storage)

##########################################################################
    # Optimise the energy system and plot the results
    ##########################################################################

    logging.info("Optimise the energy system")

    # initialise the operational model
    model = Model(energysystem)
    
    # This is for debugging only. It is not(!) necessary to solve the problem
    # and should be set to False to save time and disc space in normal use. For
    # debugging the timesteps should be set to 3, to increase the readability
    # of the lp-file.
    if debug:
        filename = os.path.join(
            helpers.extend_basic_path("lp_files"), "basic_example.lp"
        )
        logging.info("Store lp-file in {0}.".format(filename))
        model.write(filename, io_options={"symbolic_solver_labels": False})

    # if tee_switch is true solver messages will be displayed
    logging.info("Solve the optimization problem")
    model.solve(solver=solver, solve_kwargs={"tee": solver_verbose})

    logging.info("Store the energy system with the results.")

    # The processing module of the outputlib can be used to extract the results
    # from the model transfer them into a homogeneous structured dictionary.
    
    # add results to the energy system to make it possible to store them.
    energysystem.results["main"] = processing.results(model)
    energysystem.results["meta"] = processing.meta_results(model)

    # The default path is the '.oemof' folder in your $HOME directory.
    # The default filename is 'es_dump.oemof'.
    # You can omit the attributes (as None is the default value) for testing
    # cases. You should use unique names/folders for valuable results to avoid
    # overwriting.

    # store energy system with results
    energysystem.dump(dpath=None, filename=None)
    
    # *************************************************************************
    # ********** PART 2 - Processing the results ******************************
    # *************************************************************************

    logging.info("**** The script can be divided into two parts here.")
    logging.info("Restore the energy system and the results.")
    energysystem = EnergySystem()
    energysystem.restore(dpath=None, filename=None)

    # define an alias for shorter calls below (optional)
    results = energysystem.results["main"]
    storage = energysystem.groups["storage"]

    # print a time slice of the state of charge
    print("")
    print("********* State of Charge (slice) *********")
    print(
        results[(storage, None)]["sequences"][
            datetime(2024, 1, 1, 0, 0, 0) : datetime(2024, 1, 8, 0, 0, 0)
        ]
    )
    print("")

    # get all variables of a specific component/bus
    custom_storage = views.node(results, "storage")
    electricity_bus = views.node(results, "electricity")

    # plot the time series (sequences) of a specific component/bus

    fig, ax = plt.subplots(figsize=(10, 5))
    custom_storage["sequences"].plot(
        ax=ax, kind="line", drawstyle="steps-post"
    )
    plt.legend(
        loc="upper center",
        prop={"size": 8},
        bbox_to_anchor=(0.5, 1.25),
        ncol=2,
    )
    fig.subplots_adjust(top=0.8)
    plt.show()

    fig, ax = plt.subplots(figsize=(10, 5))
    electricity_bus["sequences"].plot(
        ax=ax, kind="line", drawstyle="steps-post"
    )
    plt.legend(
        loc="upper center", prop={"size": 8}, bbox_to_anchor=(0.5, 1.3), ncol=2
    )
    fig.subplots_adjust(top=0.8)
    plt.show()

   # Convert the flow results to DataFrames
    df_custom_storage = pd.DataFrame(custom_storage["sequences"])
    df_electricity_bus = pd.DataFrame(electricity_bus["sequences"])

    # Save the DataFrames to CSV files

df_custom_storage.to_csv("custom_storage_flows.csv")

df_electricity_bus.to_csv("electricity_bus_flows.csv")

    # print the solver results
    print("********* Meta results *********")
    pp.pprint(energysystem.results["meta"])
    print("")

    # print the sums of the flows around the electricity bus
    print("********* Main results *********")
    print(electricity_bus["sequences"].sum(axis=0))

    # Ergebnisse in csv-Speichern
    
    # Metaergebnisse in DataFrame konvertieren
    meta_results = energysystem.results["meta"]
    df_meta = pd.DataFrame(meta_results.items(), columns=["Parameter", "Value"])
    
    # CSV-Datei für Metaergebnisse speichern
    df_meta.to_csv("meta_results.csv", index=False)
    
    # Hauptergebnisse in DataFrame konvertieren
    main_results = electricity_bus["sequences"].sum(axis=0)
    df_main = pd.DataFrame({"Time": main_results.index, "Sum of Flows": main_results.values})
    
    # CSV-Datei für Hauptergebnisse speichern
    df_main.to_csv("main_results.csv", index=False)

if __name__ == "__main__":
    main()

I would greatly appreciate any help, particularly with implementing the cost function and integrating the solar storage system. If anyone has experience with similar issues or suggestions for solving these tasks, I would appreciate any feedback!

Thank you in advance for your assistance.

Best regards,
Armend
CSV.zip (3.0 KB)

Hi Armend,

  1. I guess the cost function that you are looking for is already implemented in oemof. The cost function is defined based on all the costs of the specified components (in your case only variable costs). This function then is to be minimized by the solver. You afterwards can a) read the result (=objective) from the solver object to determine your overall costs or b) retrieve all relevant flows from the solved model and multiply them with the corresponding cost vector. (I would suggest b) for more transparency)
  2. Simply deactivate your PV system and run your model again. Additionally I suggest a PV system module consisting of a bus with a source (PV array) and sink (excess) which then is connected to the main system (via a converter or directly). The main system could consist of a DC bus with the storage which then is connected to an AC bus using two converters (AC/DC and DC/AC).
  3. Use fixed variable costs for the supplier instead of the csv with flexible electricity prices.
    Furthermore: It will help to avoid spaghetti code and structure your code into different function or even choose an object oriented approach.

Best regards
Brian

1 Like

I have completed point 2 so far, is that correct?
Is there an example for point 1b on how to do this best?
Are these the total costs? (Image attached)
cost
CSV_DATA.zip (3.2 KB)


###########################################################################
# imports
###########################################################################
import pandas as pd
import matplotlib.pyplot as plt

import logging
import os
import pprint as pp

from datetime import datetime
from oemof.tools import logger
from oemof.solph import flows

from oemof.solph import EnergySystem
from oemof.solph import Model
from oemof.solph import buses
from oemof.solph import components as cmp

from oemof.solph import helpers
from oemof.solph import processing
from oemof.solph import views

from oemof.solph.components.experimental import SinkDSM # Corrected import statement

def main():
    # *************************************************************************
    # ********** PART 1 - Define and optimise the energy system ***************
    # *************************************************************************

    # Read data file

    verbrauchsdaten_mit_pv_file = os.path.join(os.getcwd(), 'verbrauchsdaten_mit_pv.csv')
    strom_co2_preise_file = os.path.join(os.getcwd(), 'strom_co2_preise_täglich.csv')

    verbrauchsdaten_mit_pv_data = pd.read_csv(verbrauchsdaten_mit_pv_file)
    strom_co2_preise_data = pd.read_csv(strom_co2_preise_file)

    solver = "cbc"  # 'glpk', 'gurobi',....
    debug = False  # Set number_of_timesteps to 3 to get a readable lp-file.
    solver_verbose = False  # show/hide solver output

    # initiate the logger (see the API docs for more information)
    logger.define_logging(
        logfile="oemof_example.log",
        screen_level=logging.INFO,
        file_level=logging.INFO,
    )

    date_time_index = pd.date_range('1/1/2024 00:00:00', periods=48, freq='30min')

    # create the energysystem and assign the time index
    energysystem = EnergySystem(
        timeindex=date_time_index, infer_last_interval=False
    )

    ##########################################################################
    # Create oemof objects
    ##########################################################################

    logging.info("Create oemof objects")

    # The bus objects were assigned to variables which makes it easier to
    # connect components to these buses (see below).

    # create electricity bus (AC)
    bel = buses.Bus(label="electricity")

    # create DC bus
    b_dc = buses.Bus(label="dc_bus")

    # create PV bus
    b_pv = buses.Bus(label="pv_bus")

    # adding the buses to the energy system
    energysystem.add(bel, b_dc, b_pv)

    # create excess component for the electricity bus to allow overproduction
    excess_bel = cmp.Sink(
        label="excess_bel",
        inputs={bel: flows.Flow(variable_costs=1 * -1)},  # strom verkauf = *-1
    )
    energysystem.add(excess_bel)

    # create excess component for the PV bus to allow overproduction
    excess_pv = cmp.Sink(
        label="excess_pv",
        inputs={b_pv: flows.Flow(variable_costs=1)},  # strom verkauf = *-1
    )
    energysystem.add(excess_pv)

    # create source object representing the supplier
    supplier = cmp.Source(
        label="supplier",
        outputs={
            bel: flows.Flow(
                nominal_value=1,  # 1.000 W = hier 1kwh
                variable_costs=5, # Hier sind die variablen Kosten in Euro pro Einheit
                custom_attributes={"emission_factor": 0.20},
            )
        },
    )
    energysystem.add(supplier)
    
    # create source object representing pv power plants
    pv = cmp.Source(
        label="pv",
        outputs={b_pv: flows.Flow(fix=verbrauchsdaten_mit_pv_data["pv"], nominal_value=1)},
    )  # werte in csv sind in kW also x 1000 = W
    energysystem.add(pv)
    
    # create sink object representing the electrical demand 1
    consumer1 = cmp.Sink(
        label="consumer1",
        inputs={
            bel: flows.Flow(
                fix=verbrauchsdaten_mit_pv_data["consumer1"],
                nominal_value=0.001,
                custom_attributes={"emission_factor": 0.60},
            )
        },
    )  # bedeutet wert in csv x 1 W # Beispiel: Variable Kosten von 0,25 €/kWh
    energysystem.add(consumer1)
    
    # create sink object representing the electrical demand 2
    consumer2 = cmp.Sink(
        label="consumer2",
        inputs={
            bel: flows.Flow(
                fix=verbrauchsdaten_mit_pv_data["consumer2"],
                nominal_value=0.001,
                custom_attributes={"emission_factor": 0.30},
            )
        },
    )  # bedeutet wert in csv x 1 # Beispiel: Variable Kosten von 0,25 €/kWh
    energysystem.add(consumer2)

    # create storage object representing a battery
    storage = cmp.GenericStorage(
        nominal_storage_capacity=25,# 25.000 Wh = 25 kWh oder 10,08 MWh
        label="storage",
        inputs={b_dc: flows.Flow(nominal_value=25 / 6)},# 10,08 Mwh/6 = 1,68 MW
        outputs={
            b_dc: flows.Flow(nominal_value=25 / 6, variable_costs=0.001) # kosten von 0.001 pro entladen Wh
        },
        loss_rate=0.00,
        initial_storage_level=0,
        balanced=True,
        inflow_conversion_factor=1,
        outflow_conversion_factor=0.8,
    )

    energysystem.add(storage)

    # create AC/DC converter
    ac_dc_converter = cmp.Converter(
        label="ac_dc_converter",
        inputs={bel: flows.Flow()},
        outputs={b_dc: flows.Flow()},
        conversion_factors={b_dc: 0.95}  # 95% efficiency
    )
    energysystem.add(ac_dc_converter)

    # create DC/AC converter
    dc_ac_converter = cmp.Converter(
        label="dc_ac_converter",
        inputs={b_dc: flows.Flow()},
        outputs={bel: flows.Flow()},
        conversion_factors={bel: 0.95}  # 95% efficiency
    )
    energysystem.add(dc_ac_converter)

    # connect PV bus to DC bus via a converter (assuming 100% efficiency for simplicity)
    pv_dc_converter = cmp.Converter(
        label="pv_dc_converter",
        inputs={b_pv: flows.Flow()},
        outputs={b_dc: flows.Flow()},
        conversion_factors={b_dc: 1.0}  # 100% efficiency
    )
    energysystem.add(pv_dc_converter)






    ##########################################################################
    # Optimise the energy system and plot the results
    ##########################################################################

    logging.info("Optimise the energy system")

    # initialise the operational model
    model = Model(energysystem)

    # This is for debugging only. It is not(!) necessary to solve the problem
    # and should be set to False to save time and disc space in normal use. For
    # debugging the timesteps should be set to 3, to increase the readability
    # of the lp-file.
    if debug:
        filename = os.path.join(
            helpers.extend_basic_path("lp_files"), "basic_example.lp"
        )
        logging.info("Store lp-file in {0}.".format(filename))
        model.write(filename, io_options={"symbolic_solver_labels": False})

    # if tee_switch is true solver messages will be displayed
    logging.info("Solve the optimization problem")
    model.solve(solver=solver, solve_kwargs={"tee": solver_verbose})

    logging.info("Store the energy system with the results.")

    # The processing module of the outputlib can be used to extract the results
    # from the model transfer them into a homogeneous structured dictionary.

    # add results to the energy system to make it possible to store them.
    energysystem.results["main"] = processing.results(model)
    energysystem.results["meta"] = processing.meta_results(model)

    # The default path is the '.oemof' folder in your $HOME directory.
    # The default filename is 'es_dump.oemof'.
    # You can omit the attributes (as None is the default value) for testing
    # cases. You should use unique names/folders for valuable results to avoid
    # overwriting.

    # store energy system with results
    energysystem.dump(dpath=None, filename=None)

    # *************************************************************************
    # ********** PART 2 - Processing the results ******************************
    # *************************************************************************

    logging.info("**** The script can be divided into two parts here.")
    logging.info("Restore the energy system and the results.")
    energysystem = EnergySystem()
    energysystem.restore(dpath=None, filename=None)

    # define an alias for shorter calls below (optional)
    results = energysystem.results["main"]
    storage = energysystem.groups["storage"]

    # print a time slice of the state of charge
    print("")
    print("********* State of Charge (slice) *********")
    print(
        results[(storage, None)]["sequences"][
            datetime(2024, 1, 1, 0, 0, 0) : datetime(2024, 1, 8, 0, 0, 0)
        ]
    )
    print("")

    # get all variables of a specific component/bus
    custom_storage = views.node(results, "storage")
    electricity_bus = views.node(results, "electricity")
    pv_bus = views.node(results, "pv_bus")
    dc_bus = views.node(results, "dc_bus")




    # plot the time series (sequences) of a specific component/bus

    fig, ax = plt.subplots(figsize=(10, 5))
    custom_storage["sequences"].plot(
        ax=ax, kind="line", drawstyle="steps-post"
    )
    plt.legend(
        loc="upper center",
        prop={"size": 8},
        bbox_to_anchor=(0.5, 1.25),
        ncol=2,
    )
    fig.subplots_adjust(top=0.8)
    plt.show()

    fig, ax = plt.subplots(figsize=(10, 5))
    electricity_bus["sequences"].plot(
        ax=ax, kind="line", drawstyle="steps-post"
    )
    plt.legend(
        loc="upper center", prop={"size": 8}, bbox_to_anchor=(0.5, 1.3), ncol=2
    )
    fig.subplots_adjust(top=0.8)
    plt.show()

    # electricity_bus_flows

    # Convert the flow results to a DataFrame
    df_electricity_bus = pd.DataFrame(electricity_bus["sequences"])

    # Reset the index and rename the column
    df_electricity_bus.reset_index(inplace=True)
    df_electricity_bus.rename(columns={"index": "Timestep"}, inplace=True)

    # Save the DataFrame to a CSV file
    df_electricity_bus.to_csv("electricity_bus_flows.csv", index=False)

    # custom_storage_flows

    # Convert the flow results to a DataFrame
    df_custom_storage = pd.DataFrame(custom_storage["sequences"])

    # Save the DataFrame to a CSV file with the first column named "Timestep"
    df_custom_storage.to_csv("custom_storage_flows.csv", index_label="Timestep")

    # print the solver results
    print("********* Meta results *********")
    pp.pprint(energysystem.results["meta"])
    print("")

    # print the sums of the flows around the electricity bus
    print("********* Main results *********")
    print(electricity_bus["sequences"].sum(axis=0))

    # Ergebnisse in csv-Speichern

    # Metaergebnisse in DataFrame konvertieren
    meta_results = energysystem.results["meta"]
    df_meta = pd.DataFrame(meta_results.items(), columns=["Parameter", "Value"])

    # CSV-Datei für Metaergebnisse speichern
    df_meta.to_csv("meta_results.csv", index=False)

    # Hauptergebnisse in DataFrame konvertieren
    main_results = electricity_bus["sequences"].sum(axis=0)
    df_main = pd.DataFrame({"Time": main_results.index, "Sum of Flows": main_results.values})

    # CSV-Datei für Hauptergebnisse speichern
    df_main.to_csv("main_results.csv", index=False)
    
    

if __name__ == "__main__":
    main()

Hi @Armend

To my understanding the displayed costs are all costs caused by variable costs or investments during your simulaton period.

Regarding the optimized flows:

results = solph.processing.results(model)  # get results from model
flow = results[(outflow, ac_bus)]['sequences']['flow']  # use your components instead
size = results[(outflow, ac_bus)]['scalars']['invest']  # use your components instead

Hope this helps

Regards
Brian

1 Like