PEM Electrolyzer Modeling

Q&A
, , ,
1

Hi everyone!
I am working on a problem to model a hydrogen system with a PEM electrolyzer, storage and fuel synthesis unit in Oemof. I want to model the part-load behavior of an electrolyzer using an offset transformer or a Piece-wise linear Transformer? My goal is also to determine the optimal capacity using the investment mode. Which kind of Transformer is better suited for this application?
Currently, I am using piece-wise linear transformer to model this electrolyzer. I am confused on how to model the in-break points and conversion factor so that for input electricity at each corresponding load point and efficiency, the output hydrogen is produced.
After defining sources, sinks, max. capacity for electrolyzer (testing without investment mode), is the following approach correct? (I found it in SMOOTH, that is based on oemof):


# Define the maximum capacity of the electrolyzer
max_capacity = 100e6  # 100 MW
 # Define the load points and corresponding efficiencies
bp_load_h2_prod = [0, 0.10, 0.20, 0.30, 0.40, 0.50, 0.60, 0.70, 0.80, 0.90, 1]
bp_eff_h2_prod = [0, 0.10, 0.50, 0.66, 0.71, 0.67, 0.63, 0.585, 0.545, 0.52, 0.51]
lower_heating_value = 33.3  # lower heating value in kWh/kg

bp_elec_consumed_h2_prod = [this_bp * max_capacity
                            for this_bp in bp_load_h2_prod]

bp_h2_production = []
for i_bp in range(len(bp_load_h2_prod)):
    this_hydrogen_production = \
        (bp_elec_consumed_h2_prod[i_bp]
         * bp_eff_h2_prod[i_bp]) / (lower_heating_value)
    bp_h2_production.append(this_hydrogen_production)

def get_h2_production_by_electricity(electricity_consumption):
        """Gets the hydrogen produced by the according electricity consumption
        value.
        :param electricity_consumption: electrcity consumption value [kWh]
        :return: according hydrogen production value [kg]
        """
        # Check the index of this load point.
        this_index = bp_elec_consumed_h2_prod.index(electricity_consumption)
        return bp_h2_production[this_index]

# Create and add PiecewiseLinearTransformer

pem_ely = solph.components.experimental.PiecewiseLinearTransformer(
    label="pem_ely",
    inputs={b_el: solph.flows.Flow(nominal_value=100000, variable_costs=1)},
    outputs={b_h2: Flow()},
    in_breakpoints=bp_elec_consumed_h2_prod,
    conversion_function=get_h2_production_by_electricity,
    pw_repn='CC')
energysystem.add(pem_ely)

@pschoen I see that you are currently active on oemof posts. May be you can comment on this!

Thanks in advance!
Nouman

2

Hi @nomi,

Thanks for reaching out. Honestly, I find it hard to read your code example, as you did not format it as code. Also, I am not the most experienced person on mixed-integer formulations. Thus, I’ll stay at a rather abstract level: As an experimental component, the PiecewiseLinearTransformer is poorly documented. So, if the OffsetTransformer is sufficiently detailed, I’d recommend to use that one. However, currently the size of both cannot be optimised. There is Pull Request #895, Upgrade OffsetConverter by sasa821 to change that. So, either you use this in-development feature or you will have to try different sizes until you find the optimal one.

Cheers,
Patrik

3

@nomi: on the question of code formatting, you can surround your code block in triple backslashes and optionally add the word python (or whichever language you are using) following the first set of backslashes to make your information more readable. But you still have to indent your code correctly for yourself. More here:

And the relevant example is here (although adding javascript would normally be better practice):

It is considered good netiquette to make your code snippets clean and readable (as @pschoen indicated). Also, you are more likely to get sensible answers if people can more easily see and understand what you are trying to convey. R

4

Hi @robbie.morrison Thanks for the input. I have corrected my code snippet and made it more readable.

@pschoen Thanks for your input. I hope the code is more readable now.

The concept behind the code is similar to that described here: Link

I want to implement the following curve for electrolyzer:

part_load_ely
part_load_ely1179×765 13.9 KB

Do you have any other input with edited code and this curve?

I will give a try to in-development feature for offset transformer to implement this part-load efficiency curve and will comment further on this matter.

best regards,
Nouman

1 Like
5

The code you posted seems to implement the desired curve for a constant size. I think, SMOOTH does the sizing, as this kind of optimisation is not implemented yet for the PiecewiseLinearConverter.

It might not the most elegant or efficient way to model this, but I see the following approach to solve you problem using only existing functionality of solph:

  • Add several OffsetConverters with invest optimisation (from the upcoming v0.5.1 feature that is already implemented in the aforementioned PR). If you are willing to sacrifice accuracy for performance, you could e.g. have one slope from 0 % to 20 % part load, one for 20 % to 35 % and one from 35 % to 100 %.
  • Use oemof.solph.constraints.equate_variables to set all the optimised sizes to be equal.
  • Use oemof.solph.constraints.limit_active_flow_count to make sure just one of the OffsetConverters is active.

PS: Please note that I already use the new naming scheme, as Transformer will be renamed to Converter soon.

1 Like
6

@pschoen when will v0.5.1 be out? Currently, I am not able to use OffsetTransformers or Offset Converters with investment optimzation. Is this correct or not?

BR
Nouman

7

Yes, that feature is planned for v0.5.1. It is implemented and almost ready to be included. For the release date, honestly, it’s hard to tell. We are past the date we originally planned. So, probably everything not mature will be delayed to the next release and we will have a short release sprint soon.

1 Like
8

Hi @pschoen,
I implemented the approach with several OffsetConverters(for representation of part-load behaviour of PEM Electrolyzer) with investment optimization using the dev. version.
I implemented 3 OffsetConverters and also used both of the constraints you suggested.
But I am getting the following error:

ERROR: Constructing component ‘active_electrolyzer_limit_constraint_build’
from data=None failed: KeyError: ‘Index '(“<oemof.solph.buses._bus.Bus:
'electricity_0'>”, “<oemof.solph.buses._bus.Bus: 'hydrogen_0'>”, 0)' is
not valid for indexed component 'InvestNonConvexFlowBlock.status'’

How to resolve this error or approach this problem?

Here is a code snippet for one of the OffsetConverters:

# offset_converter_2 with a minimum load of 20% and a maximum load of 35%

min_load_1 = 0.2
max_load_1 = 0.35
min_efficiency_1 = 0.50
max_efficiency_1 = 0.70

c1_1 = (max_load_1 / max_efficiency_1 - min_load_1 / min_efficiency_1) / (
        max_load_1 - min_load_1)
c0_1 = min_load_1 * (1 / min_efficiency_1 - c1_1)

pem_electrolyzer_1 = solph.components.OffsetConverter(
    label="pem_electrolyzer_1",
    inputs={b_el_1: solph.flows.Flow()},
    outputs={b_h2_1: solph.flows.Flow(nominal_value=None,
                                      variable_costs=variable_cost_pem,
                                      min=min_load_1,
                                      max=max_load_1,
                                      nonconvex=solph.NonConvex(),
                                      investment=solph.Investment(
                                          ep_costs=epc_electrolyzer,
                                          maximum=100000,
                                      ))},
    coefficients=(c0_1, c1_1),
)

Following is the code for equate_variables and limit_active_flow_count:

pem0 = energysystem.groups['pem_electrolyzer_0']
pem1 = energysystem.groups['pem_electrolyzer_1']
pem2 = energysystem.groups['pem_electrolyzer_2']
solph.constraints.equate_variables(
    optimization_model,
    optimization_model.InvestNonConvexFlowBlock.invest[pem0, b_h2_0, 0],
    optimization_model.InvestNonConvexFlowBlock.invest[pem1, b_h2_1, 0],
    optimization_model.InvestNonConvexFlowBlock.invest[pem2, b_h2_2, 0])

# limit_active_flow_count constraint
optimization_model = limit_active_flow_count(
    model=optimization_model,
    constraint_name="active_electrolyzer_limit",
    flows=[(b_el_0, b_h2_0), (b_el_1, b_h2_1), (b_el_2, b_h2_2)],
    lower_limit=0,
    upper_limit=1
)

Best regards,
Nouman

9

Hi again,

I think the problem is in the following part:

The Flows (b_el_{n}, b_h2_{n}) do not exist, instead there are (b_el_{n}, pem_electrolyzer_{n}) and (pem_electrolyzer_{n}, b_h2_{n}). As the latter is the non-convex one, you need to use it in the constraint.

10

Hi,

Thanks @pschoen for identifying the problem. That problem is now resolved. I have implemented the offset converters and the constraints to a test code. The constraints are now working to an extent.

Under variable electricity price input and H2 demand, the model always seems to consider using the 3rd offset converter for the output and does not consider storage investment in most of the cases. Although, it tends to use the storage when the load is less instead of using the other two offset converters. However, for my main use case where the fuel production is constant throughout the year, and the optimizer should choose between the times when the electricity price is lower or higher, but it tends to operate at the same load even during times with peak electricity price without using a storage.

The other problem seems to be the electricity input (which is zero sometimes for the third offset converter to my surprise even though there is an output generated) which is confusing me since it should be dependent on the load and the efficiency for the respective load point.

Can anyone shed a light on this on where I might be going wrong with this approach.

Probably, it has something to do with the conversion factor. When I initially implemented PEM as a converter with fixed efficiency, i used the following conversion factor; b_h2: efficiency/33.33, (for conversion from kW to kg of hydrogen). I am not sure how I should implement this for an offset converter.

I am using following offsetconverters:


# -------------------- Converters --------------------

# Multiple offset converters for representing piecewise linear function

# offset_converter_1 with a minimum load of 0 and a maximum load of 0.2

min_load_1 = 0.0
max_load_1 = 0.2
min_efficiency_1 = 0.0
max_efficiency_1 = 0.5
c1_1 = 2.5
c0_1 = 0

pem_electrolyzer_1 = solph.components.OffsetConverter(
    label="pem_electrolyzer_1",
    inputs={b_el_1: solph.flows.Flow()},
    outputs={b_h2_1: solph.flows.Flow(nominal_value=None,
                                      variable_costs=variable_cost_pem,
                                      min=min_load_1,
                                      max=max_load_1,
                                      investment=solph.Investment(
                                          ep_costs=epc_electrolyzer,
                                          maximum=1000,
                                      ),
                                      nonconvex=solph.NonConvex(),
                                      )},
    coefficients=(c0_1, c1_1),
)

# offset_converter_2 with a minimum load of 0.2 and a maximum load of 0.35

min_load_2 = 0.2
max_load_2 = 0.4
min_efficiency_2 = 0.50
max_efficiency_2 = 0.71
c1_2 = 1.05
c0_2 = 0.29

pem_electrolyzer_2 = solph.components.OffsetConverter(
    label="pem_electrolyzer_2",
    inputs={b_el_2: solph.flows.Flow()},
    outputs={b_h2_2: solph.flows.Flow(nominal_value=None,
                                      variable_costs=variable_cost_pem,
                                      # full_load_time_max=8000,
                                      # full_load_time_min=4000,
                                      min=min_load_2,
                                      max=max_load_2,
                                      investment=solph.Investment(
                                          ep_costs=epc_electrolyzer,
                                          maximum=1000,
                                      ),
                                      nonconvex=solph.NonConvex(),
                                      )},
    coefficients=(c0_2, c1_2),
)

# offset_converter_3 with a minimum load of 0.35 and a maximum load of 1
min_load_3 = 0.4
max_load_3 = 1
min_efficiency_3 = 0.71
max_efficiency_3 = 0.51
c1_3 = -0.3333
c0_3 = 0.8433

pem_electrolyzer_3 = solph.components.OffsetConverter(
    label="pem_electrolyzer_3",
    inputs={b_el_3: solph.flows.Flow()},
    outputs={b_h2_3: solph.flows.Flow(nominal_value=None,
                                      variable_costs=variable_cost_pem,
                                      min=min_load_3,
                                      max=max_load_3,
                                      investment=solph.Investment(
                                          ep_costs=epc_electrolyzer,
                                          maximum=1000,
                                      ),
                                      nonconvex=solph.NonConvex(),
                                      )},
    coefficients=(c0_3, c1_3),
)

# Collecting hydrogen from individual buses to the main hydrogen bus

for bus in [b_h2_1, b_h2_2, b_h2_3, b_h2_sh]:
    energysystem.add(solph.components.Converter(label=f"collect_{bus.label}",
                                                inputs={bus: solph.Flow()},
                                                outputs={b_h2: solph.Flow()},
                                                conversion_factors={bus: 1}))

and here the constraints again

optimization_model = solph.Model(energysystem)

pem1 = energysystem.groups['pem_electrolyzer_1']
pem2 = energysystem.groups['pem_electrolyzer_2']
pem3 = energysystem.groups['pem_electrolyzer_3']
solph.constraints.equate_variables(
    optimization_model,
    optimization_model.InvestNonConvexFlowBlock.invest[pem1, b_h2_1, 0],
    optimization_model.InvestNonConvexFlowBlock.invest[pem2, b_h2_2, 0],
    optimization_model.InvestNonConvexFlowBlock.invest[pem3, b_h2_3, 0])

# Apply the constraint to your model
optimization_model = limit_active_flow_count(
    model=optimization_model,
    constraint_name="active_electrolyzer_limit",
    flows=[(pem_electrolyzer_1, b_h2_1), (pem_electrolyzer_2, b_h2_2), (pem_electrolyzer_3, b_h2_3)],
    lower_limit=0,
    upper_limit=1,
)

Thanks in advance.

Best regards,
Nouman

11

First of all, the resulting efficiency curve looks reasonable, at least there are no obvious errors:

image

For the model, I do not get why you introduce the b_h2_{x} instead of just connecting the electrolyzers directly to b_h2. However, I cannot judge why you experience unexpected results from the model. It might just be that investing into the storage is not worth it. (Just observing the load curve, I would expect that also the orange part is used some times.)

1 Like
12

Thanks again for the input. The resulting efficiency curve indeed looks correct. I think the unexpected results are due to some input parameters and conversion factors taken into account. Otherwise, the model seems to be working correctly.
However, I am not sure about the role of co-efficients c1 and c0 for input-output conversion. Right now it is hydrogen output = electricity input * c1, which is problematic since it should be the efficiency at the corresponding load point which should be multiplied with the electricity input to have a suitable output.

Best regards,
Nouman

13

According to the code, it is hydrogen output = electricity input * c1 + status * c0. (Unfortunately, the documentation seems to be wrong here.) With the corresponding minimum and maximum loads, you will find the load curve as described above.