Enforce nonconvex for OffsetConverter on a single flow instead of outputs

fwitte commented 4 months ago

Resolve #1010

pep8speaks commented 4 months ago

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! :beers:

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

p-snft commented 4 months ago

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 commented 4 months ago

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 commented 4 months ago

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 commented 4 months ago

* 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 commented 4 months ago

Example 1

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

offset-two-output

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.

offset-single-output

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 commented 4 months ago

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 commented 4 months ago

I added two "convenience functions", i.e. calculate_slope_and_offset_with_reference_to_ouput and calculate_slope_and_offset_with_reference_to_input. You provide the output/input nominal power, minimal power as well as the efficiencies at both points and it will give you the slope and the offset to describe those points. I think, it might be helpful for the users to have them available within oemof-solph. On the other hand, they somehow work like a preprocessing step for the OffsetConverter, therefore they could also be integrated as an own component masquerading/facading the OffsetConverter as alternative.

What do you think, is there a place for this kind of function inside oemof-solph?

fwitte commented 4 months ago

@oemof-developer: This PR is ready for review. It breaks with the API of the last version(s). In the new version, the following things changed:

The OffsetConverter is now multi-input and multi-output, before only multi-output was possible.
The offset and the slope are now not provided with the coefficients keyword, but by conversion_factors (slope) and normed_offsets (offset)
There must only be a single NonConvex flow in the OffsetConverter. In case of investment optimization, the investment has to be provided to the same flow. Before, any one of the outputs could have the NonConvex attribute, I am not sure, how a multi-output OffsetConverter was handled with investment optimization simultaneously.
(It seems?) the nominal_value attribute has to be provided to the NonConvex flow.
All offsets and slopes reference the flow and status_nominal of the NonConvex flow. The offset always referenced the flow with status_nominal, while the slope could reference the input flow.

fwitte commented 4 months ago

Hi @Bachibouzouk,

Thanks @fwitte for this PR :) Nice work! The review I made is not very in depth, I haven't tested it locally because of a lack of example I could readily run (i.e. lack of time to build such an example). However I still have a few suggestions that might be relevant:

Thank you for the feedback and your time :)

* an example like https://github.com/oemof/oemof-solph/blob/dev/examples/offset_converter_example/offset_diesel_genset_nonconvex_investment.py would be nice

Ah, I overlooked that, thank you for the reminder. I will add one.

* what about the back compatibility? Shouldn't this component bear another name to avoid breaking changes?

How is this usually dealt with in oemof-solph? I could think of:

Add the coefficients back to the new version of the OffsetConverter and make an internal conversion into the new format.
At the same time, when users use the coefficients they will get a warning, that the old API will be deprecated with the next major release.

Or would you actually suggest giving the new component a new class name and simply phase out the OffsetConverter? In the end, the user, who wants to keep his/her source code up to date, has to update in anyway: Either by modifying the specifications for the component or by moving to the new one, or? I would be interested in your perspective.

* an amendment to the documentation of the OffsetConverter would be good

That there is now a new component, that should be used instead of the OffsetConverter?

Thank you, all the best

Francesco

Bachibouzouk commented 4 months ago

Or would you actually suggest giving the new component a new class name and simply phase out the OffsetConverter?

I made this comment in light of what @p-snft said in https://github.com/oemof/oemof-solph/issues/1010#issuecomment-2115552441

I would be interested in your perspective.

I agree with you that when component evolve the users have to adapt the code, I think for that a warning that coefficients is not used anymore would be the minimum. If from the coefficient you could compute slope and offset, then this is would be ideal. If we pursue the "create a new component" path, then it simplifies the back-compatibility issues as we just need to provide a warning on the use of OffsetConverter

That there is now a new component, that should be used instead of the OffsetConverter?

I am not aware of this, I meant more that the documentation needs to reflect the new OffsetConverter you are proposing

fwitte commented 4 months ago

Or would you actually suggest giving the new component a new class name and simply phase out the OffsetConverter?

I made this comment in light of what @p-snft said in #1010 (comment)

Thank you for the hint

I would be interested in your perspective.

I agree with you that when component evolve the users have to adapt the code, I think for that a warning that coefficients is not used anymore would be the minimum. If from the coefficient you could compute slope and offset, then this is would be ideal. If we pursue the "create a new component" path, then it simplifies the back-compatibility issues as we just need to provide a warning on the use of OffsetConverter

I'll look into it and ping to rereview once I have all other changes ready.

That there is now a new component, that should be used instead of the OffsetConverter?

I am not aware of this, I meant more that the documentation needs to reflect the new OffsetConverter you are proposing

That is actually already updated (it is the section where you asked about the l_max)

fwitte commented 3 months ago

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.

fwitte commented 3 months ago

@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 commented 3 months ago

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.

fwitte commented 3 months ago

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.

oemof / oemof-solph