Hello all,
I am working on a code to look into storage and usage of waste heat. For now, I have the following scenario which I could not find a working solution for yet:
There is a PV and/or wind power plant with a defined nominal power (2 MW) and and a load profile for a year, feeding into the electricity bus. An electrolyser with a nominal power of 1 MW is connected to said power plant, producing waste heat which is fed into the thermal bus (the produced hydrogen is not further considered here). To simplify, for now all the heat is going to a variable sink. There is also a variable sink for the electricity bus for any excess electricity not being fed into the electrolyser.
The electrolyser must only operate if the available electricity is at least 3% of its nominal power, otherwise the electrolyser must be shut down, having a thermal output of 0. This last requirement is what have not succeeded in programming yet.
What I have coded so far results in all the electricity going to the âexcess electricity sinkâ and the electrolyer not running at all. As I am rather new to optimisation my guesses are that either the constraint for the electrolyser is incorrectly set or that there is something wrong with the optimisation objective.
Code
# imports
import logging
import os
import pprint as pp
import warnings
from datetime import datetime
from datetime import time as settime
import matplotlib.pyplot as plt
import pandas as pd
from oemof.tools import logger
import numpy as np
import oemof.solph as solph
try:
import matplotlib.pyplot as plt
except ImportError:
plt = None
solver = 'cbc'
debug = False
number_of_time_steps = 8760
solver_verbose = False
# initiate the logger
logger.define_logging(
logfile="MA_logfile.log",
screen_level=logging.INFO,
file_level=logging.INFO,
)
logging.info('Initialisieren des Energiesystems')
date_time_index = pd.date_range('01/01/2018', periods=number_of_time_steps,
freq='h')
# create the energysystem and assign the time index
energysystem = solph.EnergySystem(timeindex=date_time_index, infer_last_interval=False)
###########################################################################
# Einlesen der Dateien
logging.info('Einlesen der Dateien')
filename = os.path.join(os.path.dirname(__file__), 'PV2018.csv') # Quelle: renewables ninja
PV = pd.read_csv(filename)
filename1 = os.path.join(os.path.dirname(__file__), 'wind2018.csv') # Quelle: renewables ninja
wind = pd.read_csv(filename1)
###########################################################################
# Erstellen und Hinzufuegen der Objekte zum Energiesystem
logging.info('Erstellen der Objekte')
# Wärmenetz erstellen
bth = solph.Bus(label="Waerme")
# Strombus fßr Wärmepumpe und Elektrolyseur erstellen
bel = solph.Bus(label="Strom")
# Wärme-Bus und Strom-Bus zum Energiesystem hinzufuegen
energysystem.add(bth)
energysystem.add(bel)
# Definition der Nennleistung fĂźr PV-Anlage und Windrad
nominal_power = 5000 # nominal value in kW
# PV-Anlage erstellen (fixe Source)
energysystem.add(solph.components.Source(label='PV', outputs={bel: solph.Flow(
fix=PV['electricity'], nominal_value=nominal_power)}))
# Wind-Anlage erstellen (fixe Source)
energysystem.add(solph.components.Source(label='wind', outputs={bel: solph.Flow(
fix=wind['electricity'], nominal_value=nominal_power)}))
# Definition der Nennleistung der Elektrolyseure (1 MW)
nominal_power_EL = 1000
# Definieren des Umwandlungsfaktors fßr Abwärme
conversion_factor_EL = 0.18
# Umwandlungsfaktor fßr Abwärme
conversion_factor_heat = 0.18 # 18% Umwandlungsrate von Strom in nutzbare Abwärme
# Minimum Input Threshold: 3% der Nennleistung des Elektrolyseurs
minimum_input_threshold = 0.03 * nominal_power_EL
# Elektrolyseur als Converter erstellen
electrolyser = solph.components.Converter(
label='electrolyser',
inputs={bel: solph.Flow()},
outputs={bth: solph.Flow(
nominal_value=nominal_power_EL * conversion_factor_heat,
min=0)},
conversion_factors={bth: conversion_factor_heat}
)
# Komponenten zum Energysystem hinzufĂźgen
energysystem.add(electrolyser)
# Optimierungsmodell erstellen
model = solph.Model(energysystem)
# Constraint hinzufĂźgen, das den Elektrolyseur nur aktiviert, wenn der Input Ăźber dem Minimum Threshold liegt
def minimum_input_rule(model):
for t in model.TIMESTEPS:
input_flow = model.flow[(electrolyser, solph.Bus(label='electricity')), t]
model.minimum_input_constraint = Constraint(
expr=(input_flow >= minimum_input_threshold)
)
# Constraint dem Modell hinzufĂźgen
model.minimum_input_constraint = minimum_input_rule
# Abwaerme erstellen (variable Sink)
energysystem.add(solph.components.Sink(label='Abwaerme', inputs={bth: solph.Flow()}))
# Ăberschussstrom erstellen (Sink)
energysystem.add(solph.components.Sink(label='Ăberschussstrom', inputs={bel: solph.Flow()}))
# Ergebnisse verarbeiten
results = solph.processing.results(model)
###########################################################################
# Optimierung des Energiesystems
logging.info('Optimierung des Energiesystems')
# Initialisieren des operativen Modells
model = solph.Model(energysystem)
logging.info('Loesen des Optimierungsproblems')
model.solve(solver=solver, solve_kwargs={'tee': solver_verbose})
# Hinzufuegen der Ergebnisse zum Energiesystem, um das Speichern zu ermoeglichen
energysystem.results = solph.processing.results(model)
logging.info('Wiederherstellen des Energiesystems und der Ergebnisse')
energysystem = solph.EnergySystem()
energysystem.restore(dpath=None, filename='Beispiel.oemof')
# Definition eines Kurznamens fuer die Ergebnisse
results = energysystem.results
# Variablen fuer die Ausgabe definieren
PV = solph.views.node(results, 'PV')
Abwaerme = solph.views.node(results, 'Abwaerme')
Ueberschussstrom = solph.views.node(results, 'Ăberschussstrom')
###########################################################################
# Plotten der Ergebnisse,
if plt is not None:
fig, ax = plt.subplots(figsize=(10, 5))
Abwaerme['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)
fig, ax.set_ylabel('Leistung in kW')
plt.show()
fig, ax = plt.subplots(figsize=(10, 5))
Ueberschussstrom['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)
fig, ax.set_ylabel('Leistung in kW')
plt.show()
Does anyone have ideas on how to code this?
Any help and advice is highly appreciated!
As I am working on storing and using waste heat, the temperature is also a relevant factor that I was not able to integrate into the model (using oemof.thermal) yet but I think that is a topic for a different
A municipal application.
The discussion on thermal sub-networks might be of particular interest, §8.6 (p146â150) (numbered 161 on PDF).