Model an electrolyser that has a minimum load it can operate at

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

Just noting the use of a collapsible code block in the above posting. Using markdown:

[details="Collapsible code block"]
``` python
import warnings
```
[/details]

Which renders as:

Collapsible code block
import warnings

You can also do the same thing with HTML <details> tags, see: https://meta.discourse.org/t/code-blocks-within-details

▢

1 Like

Hello @ella  On the topic of storing and using recovered heat, the following admittedly old publications might offer some insights still. Always useful to see what others have done when attempting to answer similar questions. No electrolyzers back then though.

  • Groscurth, Helmuth-M, Thomas Bruckner, and Reiner KĂźmmel (1995). “Modeling of energy-services supply systems”. Energy. 20 (9): 941–958. doi:10.1016/0360-5442(95)00067-Q. Describes the deeco framework in conceptual terms. (And which I presume is quite similar to oemof?)

  • Lindenberger, Dietmar, Thomas Bruckner, Helmuth-M Groscurth, and Reiner KĂźmmel (2000). “Optimization of solar district heating systems: seasonal storage, heat pumps, and cogeneration”. Energy. 25 (7): 591–608. doi:10.1016/S0360-5442(99)00082-1. :closed_access: A municipal application.

  • Morrison, Robbie (31 July 2000). Optimizing exergy-services supply networks for sustainability — MSc thesis. Dunedin, New Zealand: Physics Department, Otago University. doi:10.5281/zenodo.13360487. :open_access: The discussion on thermal sub-networks might be of particular interest, §8.6 (p146–150) (numbered 161 on PDF).

  • Bruckner, Thomas (2001). Benutzerhandbuch ‘‘deeco’’ â€” Version 1.0 [User handbook '‘deeco’ â€” Version 1.0] (in German). Berlin, Germany: Institut fĂźr Energietechnik (IET), Technische Universität Berlin. doi:10.5281/zenodo.5148149. 239 pages. :open_access: Lists all equations and sources.

Shut-down mode requires a mixed‑integer formulation. I will leave it to the oemof crew to have first crack at answering that question. No doubt it is either currently supported or easy to implement. Good luck, R.

1 Like

Hello,

this topic is interesting to me as we are working on similar things for our energy modelling projects. Unfortunately I can’t help with the oemof issue. The best I can do is shill our tool, which may be a good fit if optimisation plays a lesser role than control of operation.

Regarding the minimum power and waste heat of electrolysers I can report a bit from a realised 1 MW electrolyser project described in:

For more recent and accurate numbers there’s a poster in German, from a recent event:

As almost all electrolysers the described plant consists of several units, each with their own power supply and several electrolysis stacks, which are dispatched from the central control of the entire plant. These units have a high minimum power fraction (around 40%). So a low power operation could only be achieved by dispatching single units close to their minimum power. For various other practical reasons (that don’t exist in Model Land), this is results in the plant rarely if ever being operated below 25% of nominal power.

The waste heat recovery of the plant is done via two routes: One for the (relatively) high temperature (55-65 °C) waste heat of the stacks and one for various other heat sources (25-35°C) such as the power electronics, cooling the room in which the stacks are, and the hydrogen compressor (the stacks provide 10-30 bar, the storage requires 300-700 bar). The latter route is probably not possible for all other plants, especially since container solutions seem to become the standard.

I hope this helps in any way.

Hello @ella,
how you describe it I would assume that a [NonConvexFlowBlock] (oemof.solph.Flow — oemof.solph 0.6.0a2 documentation) class would do the job. With that you can define min/max values of the flow as fractions of the capacity.
Be aware that this makes your optimization nonconvex which increases the computing time.
There is a good explanation in the docs (for a diesel generator, but it’s equivalent) - - check for the header Combination of Dispatch and Investment Optimisation

Hello @georgaloa,
yes, I found the nonconvex flow but cannot piece it together correctly yet for it to work.

When I replace

Collapsible code block
# 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

with

Collapsible code block
electrolyserPV = solph.components.Converter(
    label='Elektrolyse_PV',
    inputs={bel: solph.Flow()},
    outputs={bth: solph.flows.Flow(
        min=0,
        max=nominal_power_EL * conversion_factor_EL,
        nominal_value=nominal_power_EL,
        nonconvex=solph.NonConvex()
        )},
    conversion_factors={bth: conversion_factor_EL}
)
# Komponenten zum Energysystem hinzufĂźgen
energysystem.add(electrolyserPV)

# Optimierungsmodell erstellen
model = solph.Model(energysystem)

I get the same results. I am aware that the altered code contains a nonconvex flow but no further information on the nonconvexity. From the documentation I do not understand how I can include that. The parameters that I can find in the documentation in the IDE are: initial_status, minimum_uptime, minimum_downtime, maximum_startups, maximum_shutdowns, startup_costs, shutdown_costs, activity_costs, inactivity_costs, negative_gradient_limit, positive_gradient_limit.
In the documentation for NonConvexFlow mentions a _minimum_flow_constraint() which I would like to use

When I used shutdown_costs just to try it out the run time went from less than a minute to almost ten and I got the following logging info:

WARNING-Implicitly replacing the Component attribute shutdown_costs (type=<class ‘pyomo.core.base.expression.ScalarExpression’>) on block NonConvexFlowBlock with a new Component (type=<class ‘pyomo.core.base.expression.ScalarExpression’>).
This is usually indicative of a modelling error.
To avoid this warning, use block.del_component() and block.add_component().

which I also find confusing, especially the “implicitly replacing” part. How is defining a parameter implicit or does this happen somewhere else?

Any further explanations are welcome.

Regards,
Ella

@etienne.ott, thanks for the further information provided. However, my question is really only about the coding part as I have a specific electrolyser in mind for operation in the 3-100% range. (The AEM Multicore/Nexus by Enapter, working at SIZ maybe you know about it).

Indeed, the Hydrogen Terminal Braunschweig has recently been constructed intended as test stand and reasearch electrolyser plant. I’m not personally involved, but I can put you (and anyone else interested) into contact with the group running it. Although full operation has not begun yet, so there is little insight on operation to be gleaned right now.

See below, that should actually work.

Collapsible code block
electrolyserPV = solph.components.Converter(
    label='Elektrolyse_PV',
    inputs={bel: solph.Flow()},
    outputs={bth: solph.flows.Flow(
        min=0.03,
        max=1,
        nominal_value=nominal_power_EL,
        nonconvex=solph.NonConvex()
        )},
    conversion_factors={bth: conversion_factor_EL}
)

For the nonconvex block, can you share the code which is not working?
I would recommend to also study in detail the general oemof-examples (oemof-solph/examples/start_and_shutdown_costs/startup_shutdown.py at dev ¡ oemof/oemof-solph ¡ GitHub)

Collapsible code block
    pp1 = solph.components.Source(
        label="plant_min_down_constraints",
        outputs={
            bel: solph.Flow(
                nominal_value=10,
                min=0.5,
                max=1.0,
                variable_costs=10,
                nonconvex=solph.NonConvex(
                    minimum_downtime=4, initial_status=0
                ),
            )
        },
    )

    pp2 = solph.components.Source(
        label="plant_min_up_constraints",
        outputs={
            bel: solph.Flow(
                nominal_value=10,
                min=0.5,
                max=1.0,
                variable_costs=10,
                nonconvex=solph.NonConvex(minimum_uptime=2, initial_status=1),
            )
        },
    )
1 Like

Your suggested solution works. Thanks a lot!