Enforce nonconvex for OffsetConverter on a single flow instead of outputs

[fwitte]

Resolve #1010

[pep8speaks]

Hello @fwitte! Thanks for updating this PR. We checked the lines you've touched for PEP 8 issues, and found:

There are currently no PEP 8 issues detected in this Pull Request. Cheers! 🍻

Comment last updated at 2024-06-19 07:13:06 UTC

[p-snft]

Thank you. I think this issue came up as the offset is typically negative, so there would be no output without input. However, when taking into account the losses (as you did in #1067) you will always run into the issue. So, I'd also consider the current implementation broken.

As you suggested in the chat, it might make sense to change the keyword argument, when we change the implementation anyway. I'd suggest conversion_factors for the slope, as it's the same as for the (simple) Converter, and normed_offset to emphasise the fact that it's normed.

[fwitte]

Thank you. I think this issue came up as the offset is typically negative, so there would be no output without input. However, when taking into account the losses (as you did in #1067) you will always run into the issue. So, I'd also consider the current implementation broken.

Agree!

As you suggested in the chat, it might make sense to change the keyword argument, when we change the implementation anyway. I'd suggest conversion_factors for the slope, as it's the same as for the (simple) Converter, and normed_offset to emphasise the fact that it's normed.

That sounds like a good idea. I will complete the changes in this PR then. This includes

• rename the keywords
• update the docs and examples
• fix/create tests

[fwitte]

On top of that, I suggest a couple more changes:

• Rename to MultiOutputOffsetConverter, and we could add a new component called MultiInputOffsetConverter or
• enforce a single flow of all inputs and outputs to be NonConvex. This flow would be the basis for all other conversion_factors and normed_offsets. Therefore it also provides the min and max boundaries, as well as the status_nominal and also the Invest if it should be considered.

The second version would be more simple from the user perspective (you have to define a reference flow, i.e. the NonConvex one, and all given offsets and slopes refer to that, no matter if it has multi input, multi output or both, and no matter of that NonConvex flow is an input or an output. The other version is likely more easy to implement (or maintain), but I am not really sure. What do you think?

[fwitte]

* enforce a single flow of all inputs and outputs to be `NonConvex`. This flow would be the basis for all other `conversion_factors` and `normed_offsets`. Therefore it also provides the `min` and `max` boundaries, as well as the `status_nominal` and also the `Invest` if it should be considered.

I drafted this version with my latest commit. I will readd the investment part after I have created two examples which we should discuss before putting more time into the changes.

[fwitte]

Example 1

electricity input with hydrogen (efficiency declining with load) and heat (efficiency increasing with load) output, reference is the only input flow:

from oemof import solph
import numpy as np


t_steps = 11
es = solph.EnergySystem(timeindex=solph.create_time_index(year=2023, number=t_steps - 1))


eta_h2_min = 0.6
eta_h2_max = 0.5
P_in_min = 2
P_in_max = 10

P_out_min_h2 = P_in_min * eta_h2_min
P_out_max_h2 = P_in_max * eta_h2_max

slope_h2 = (P_out_max_h2 - P_out_min_h2) / (P_in_max - P_in_min)
offset_h2 = (P_out_max_h2 - slope_h2 * P_in_max) / P_in_max


eta_heat_min = 0.25
eta_heat_max = 0.3
P_out_min_heat = P_in_min * eta_heat_min
P_out_max_heat = P_in_max * eta_heat_max

slope_heat = (P_out_max_heat - P_out_min_heat) / (P_in_max - P_in_min)
offset_heat = (P_out_max_heat - slope_heat * P_in_max) / P_in_max



b_h2 = solph.Bus("hydrogen bus")
b_heat = solph.Bus("heat bus")
b_el = solph.Bus("electricity bus")

b_electrolyzer_h2 = solph.Bus("electrolyzer virtual input for h2")
b_electrolyzer_heat = solph.Bus("electrolyzer virtual input for heat")

source_el = solph.components.Source("electricity import", outputs={b_el: solph.Flow()})
sink_heat = solph.components.Sink("heat export", inputs={b_heat: solph.Flow()})
sink_h2 = solph.components.Sink("hydrogen export", inputs={b_h2: solph.Flow(fix=[0] + np.linspace(P_out_min_h2, P_out_max_h2, t_steps - 1).tolist(), nominal_value=1)})

electrolyzer = solph.components.OffsetConverter(
    label="electrolyzer",
    inputs={b_el: solph.Flow(nominal_value=P_in_max, nonconvex=solph.NonConvex(), min=P_in_min / P_in_max)},
    outputs={
        b_heat: solph.Flow(),
        b_h2: solph.Flow()
    },
    conversion_factors={
        b_heat: slope_heat,
        b_h2: slope_h2,
    },
    normed_offsets={b_heat: offset_heat, b_h2: offset_h2}
)

es.add(
    b_el, b_heat, b_h2,
    source_el, sink_h2, sink_heat,
    electrolyzer
)

om = solph.Model(es)

om.solve("gurobi")

result = solph.views.convert_keys_to_strings(om.results())

heat = result[("electrolyzer", "heat bus")]["sequences"]["flow"]
el = result[("electricity bus", "electrolyzer")]["sequences"]["flow"]
h2 = result[("electrolyzer", "hydrogen bus")]["sequences"]["flow"]

plot_lines_and_conversion_efficiencies(
    {"flows": el, "min": P_in_min, "max": P_in_max, "commodity": "electricity input"},
    {"flows": heat, "offset": offset_heat, "slope": slope_heat, "commodity": "heat output"},
    {"flows": h2, "offset": offset_h2, "slope": slope_h2, "commodity": "h2 output"},
)

Example 2

gas input and electricity output with electricity (output commodity) as reference, might not make that much sense, but to show, that it works.

from oemof import solph
import numpy as np


t_steps = 11
es = solph.EnergySystem(timeindex=solph.create_time_index(year=2023, number=t_steps - 1))


eta_el_min = 0.3
eta_el_max = 0.4
P_in_min = 2
P_in_max = 10

P_out_min_el = P_in_min * eta_el_min
P_out_max_el = P_in_max * eta_el_max

slope_gas = (P_in_max - P_in_min) / (P_out_max_el - P_out_min_el)
offset_gas = (P_in_max - slope_gas * P_out_max_el) / P_out_max_el


b_el = solph.Bus("electricity bus")
b_gas = solph.Bus("gas bus")

source_gas = solph.components.Source("gas import", outputs={b_gas: solph.Flow()})
sink_el = solph.components.Sink("electricity export", inputs={b_el: solph.Flow(fix=[0] + np.linspace(P_out_min_el, P_out_max_el, t_steps - 1).tolist(), nominal_value=1)})

chp = solph.components.OffsetConverter(
    label="chp",
    inputs={b_gas: solph.Flow()},
    outputs={
        b_el: solph.Flow(nominal_value=P_out_max_el, nonconvex=solph.NonConvex(), min=P_out_min_el / P_out_max_el)
    },
    conversion_factors={b_gas: slope_gas},
    normed_offsets={b_gas: offset_gas}
)

es.add(
    b_gas, b_el,
    source_gas, sink_el,
    chp
)

om = solph.Model(es)

om.solve("gurobi")

result = solph.views.convert_keys_to_strings(om.results())

el = result[("chp", "electricity bus")]["sequences"]["flow"]
gas = result[("gas bus", "chp")]["sequences"]["flow"]

plot_lines_and_conversion_efficiencies(
    {"flows": el, "min": P_out_min_el, "max": P_out_max_el, "commodity": "electricity output"},
    {"flows": gas, "offset": offset_gas, "slope": slope_gas, "commodity": "gas input"},
)

Plotting code

def plot_lines_and_conversion_efficiencies(x, *args):

    fig, ax = plt.subplots(2, len(args), sharex=True)

    x_min = x["min"]
    x_max = x["max"]
    x_commodity = x["commodity"]
    x = x["flows"]

    ax = ax.transpose().flatten()

    for i, y in enumerate(args):
        offset = y["offset"]
        slope = y["slope"]
        y_commodity = y["commodity"]
        y = y["flows"]

        ax[2 * i].plot([0, x_min], [offset * x_max, offset * x_max + slope * x_min], "--", color="k")
        ax[2 * i].plot([x_min, x_max], [offset * x_max + slope * x_min, (offset + slope) * x_max], color="k")
        ax[2 * i].scatter(x, y)

        ax[2 * i].set_ylabel(y_commodity)
        ax[2 * i].grid()

        ax[2 * i + 1].scatter(x, y / x)
        ax[2 * i + 1].set_ylabel(f"{y_commodity} to {x_commodity} ratio")
        ax[2 * i + 1].grid()
        ax[2 * i + 1].set_xlabel(x_commodity)

    plt.tight_layout()
    fig.savefig("figure.png")

[fwitte]

The lp file of last failing test (TestsMultiPeriodConstraint.test_offsetconverter) cannot be reproduced identically due to the reformulation of the OffsetConverter rule. With the previous implementation we get

$$0=-\dot E_\text{therm} \cdot + 0.9 \dot E_\text{gas} - 17 y_\text{therm} \cdot 100$$

With the new implementation we have to reference both, the slope (0.9) and the offset (17) to one of both flows. That means, that either the slope or the offset has to be recalculated. I therefore suggest, hand-checking the new implementation's lp and replacing the old one.

On a side note: the offset of 17 does not make a lot of sense here, I think it might be relict from the time, when the offset was not yet normed to the nominal_status.

[fwitte]

I added two more tests, one for double inputs and double outputs with the NonConvex flow at the input, and the same with the NonConvex flow at an output.

With that, I think this PR is ready for final review.

[p-snft]

Thank you for this implementation. It really fills a gap. At the moment, there are few issues left to be fixed before merge.

[fwitte]

@p-snft, @Bachibouzouk, I think, we are finally at it. I implemented a small plotting function for the new figures. I did not reference it anywhere in the docs. Although it is helpful, it could be refined a little bit maybe, I did not want to waste too much time on it. The fail of the docs test stage seems to be a false positive.

Have a nice weekend!

[p-snft]

The doc error reports that the linked page https://shop.vde.com/en/vde-study-the-cellular-approach has problems. This is no problem to be solved here. There is little to do, yet. I will contribute that myself.

[p-snft]

In my opinion, this is ready to be merged.

There are two remaining problems found by the automated process:

• Docs: Not all linked pages are available. (Not to be solved here.)
• Coverage: Test coverage decreases because of plotting functionality. As, however, this is not an advertised feature and nobody will rely on it, I think this is in order.

[fwitte]

In my opinion, this is ready to be merged.

There are two remaining problems found by the automated process:

* Docs: Not all linked pages are available. (Not to be solved here.)

* Coverage: Test coverage decreases because of plotting functionality. As, however, this is not an advertised feature and nobody will rely on it, I think this is in order.

Thank you for the review @p-snft and @Bachibouzouk, I will merge the branch although the tests are crashing here. The reason for that is that some parts of the already existing code are incompatible with numpy. I did not touch the code since the approval, where the tests were still successful.