Building EnergySystem succeeds even if Flows point out

[fwitte]

When forgetting to add elements to the EnergySystem (e.g. because you forgot) things happen, that look not intended.

Consider the following very simple system:

• Source -> Bus "b1" -> Sink

The correct code would look like this:

# Snippet 1
from oemof import solph
from oemof.solph import Bus
from oemof.solph import EnergySystem
from oemof.solph import Flow
from oemof.solph import Model
from oemof.solph.components import Source, Sink


es = EnergySystem(timeindex=solph.create_time_index(2024, number=1))


b1 = Bus("bus1")

source = Source("source", outputs={b1: Flow()})
sink = Sink("sink", inputs={b1: Flow(nominal_value=1, fix=2)})

es.add(b1, source, sink)
model = Model(es)
model.solve("cbc")

Now, if you forgot to add the sink, the result is exactly the same. The model works and the source produces what is required by the sink even though we did not add it to the EnergySystem. This is in my opinion at least surprising, but does not need to be an actual issue.

# Snippet 2
from oemof import solph
from oemof.solph import Bus
from oemof.solph import EnergySystem
from oemof.solph import Flow
from oemof.solph import Model
from oemof.solph.components import Source, Sink


es = EnergySystem(timeindex=solph.create_time_index(2024, number=1))


b1 = Bus("bus1")

source = Source("source", outputs={b1: Flow()})
sink = Sink("sink", inputs={b1: Flow(nominal_value=1, fix=2)})

es.add(b1, source)
model = Model(es)
model.solve("cbc")

However, if you were to add more components in between, e.g. change to

• Source -> Bus "b1" -> Bus "b2" -> Sink

then forgetting to add the Bus "b2" and the sink will add the flow of the Bus "b2" but will ignore the Flow attached to the sink, which seems somehow inconsistent to the behavior above. The example below would lead to an infeasibility if all flows were considered:

# Snippet 3
from oemof import solph
from oemof.solph import Bus
from oemof.solph import EnergySystem
from oemof.solph import Flow
from oemof.solph import Model
from oemof.solph.components import Source, Sink


es = EnergySystem(timeindex=solph.create_time_index(2024, number=1))


b1 = Bus("bus1")
b2 = Bus("bus2", inputs={b1: Flow(nominal_value=1, fix=1)})

source = Source("source", outputs={b1: Flow()})
sink = Sink("sink", inputs={b2: Flow(nominal_value=1, fix=2)})

es.add(b1, source)
model = Model(es)
model.solve("cbc")

Changing this to the correct version:

# Snippet 4
from oemof import solph
from oemof.solph import Bus
from oemof.solph import EnergySystem
from oemof.solph import Flow
from oemof.solph import Model
from oemof.solph.components import Source, Sink


es = EnergySystem(timeindex=solph.create_time_index(2024, number=1))


b1 = Bus("bus1")
b2 = Bus("bus2", inputs={b1: Flow(nominal_value=1, fix=1)})

source = Source("source", outputs={b1: Flow()})
sink = Sink("sink", inputs={b2: Flow(nominal_value=1, fix=2)})

es.add(b1, source, b2)  # sink again optional in the current implementation, the flow into the sink will still be considered
model = Model(es)
model.solve("cbc")

Finally, to make things even more weird: If you put the flows as inputs and outputs of the Bus "b2" (and no `inputs´ for the sink), this model works fine, even though the input and output flows on Bus "b2" contradict in their constraints:

# Snippet 5
from oemof import solph
from oemof.solph import Bus
from oemof.solph import EnergySystem
from oemof.solph import Flow
from oemof.solph import Model
from oemof.solph.components import Source, Sink


es = EnergySystem(timeindex=solph.create_time_index(2024, number=1))


b1 = Bus("bus1")

source = Source("source", outputs={b1: Flow()})
sink = Sink("sink")
b2 = Bus("bus2", inputs={b1: Flow(nominal_value=1, fix=1)}, outputs={sink: Flow(nominal_value=1, fix=2)})

es.add(b1, source)
model = Model(es)
model.solve("cbc")

The infeasibility is only detected if the Bus "b2" is actually added to the EnergySystem:

# Snippet 6
from oemof import solph
from oemof.solph import Bus
from oemof.solph import EnergySystem
from oemof.solph import Flow
from oemof.solph import Model
from oemof.solph.components import Source, Sink


es = EnergySystem(timeindex=solph.create_time_index(2024, number=1))


b1 = Bus("bus1")

source = Source("source", outputs={b1: Flow()})
sink = Sink("sink")
b2 = Bus("bus2", inputs={b1: Flow(nominal_value=1, fix=1)}, outputs={sink: Flow(nominal_value=1, fix=2)})

es.add(b1, b2, source)
model = Model(es)
model.solve("cbc")

In my opinion, it would be better, if the EnergySystem would always ignore parts of the system if not explicitly added. What do you think?

Edit by @p-snft: I allowed myself to enumerate the code snippets.

[p-snft]

I made sense out of this:

• Not adding something will not do anything (esp. not create constraints). This is expected.
• A Node has to be added explicitly. This is expected.
• A Flow is added if at least one Node it connects to is added to the energy system. This is expected.
  • Building the model will fail if the Node it leaves is not part of the EnergySystem. This is expected.
  • Building the model will still succeed if the Node it enters is not part of the EnergySystem. This is unexpected.

I think the solution would be to fail building the model if the Node a Flow points to is not part of the EnergySystem. For the implementation of the fix, I think we should have a graph structure consistency check at the level of oemof.network.

Here is a shorter version of @fwitte's test cases including the explanation what actually happens in the cases he presented. For completeness, I added two cases (0 and 7).

from oemof import solph
from oemof.solph import Bus
from oemof.solph import EnergySystem
from oemof.solph import Flow
from oemof.solph import Model

def solve_case(variation):
    es = EnergySystem(
        timeindex=solph.create_time_index(2024, number=3),
        infer_last_interval=False,
    )

    source = Bus("source", balanced=False)
    b1 = Bus("b1")
    sink = Bus("sink", balanced=False)

    b1.inputs[source] = Flow(nominal_value=2, fix=1, variable_costs=0.1)

    match variation:
        case 0:
            es.add(source, b1, sink)
            # infeasible: b1 cannot be balanced, as there is no Flow
        case 1:
            sink.inputs[b1] = Flow()  # considered via b1 and sink
            es.add(source, b1, sink)
            # optimal
        case 2:
            sink.inputs[b1] = Flow()  # still considered via b1
            es.add(source, b1)
            # optimal
        case 3:
            b2 = Bus("b2")  # not considered
            b2.inputs[b1] = Flow()  # still considered via b1
            sink.inputs[b2] = Flow(nominal_value=1)  # not considered
            es.add(source, b1)
            # optimal
        case 4:
            b2 = Bus("b2")  # considered
            b2.inputs[b1] = Flow()  # not considered
            sink.inputs[b2] = Flow(nominal_value=1)  # not considered
            es.add(source, b1, b2)
            # infeasible: b2 cannot be balanced (no Flow out)
        case 5:
            b2 = Bus("b2")  # not considered
            b2.inputs[b1] = Flow()  # considered via b1
            b2.outputs[sink] = Flow(nominal_value=1)  # not considered
            es.add(source, b1)
            # optimal
        case 6:
            b2 = Bus("b2")  # considered
            b2.inputs[b1] = Flow()  # considered
            b2.outputs[sink] = Flow(nominal_value=1)  # considered
            es.add(source, b1, b2)
            # infeasible: b2 cannot be balanced (Flow out too small)
        case 7:
            sink.inputs[b1] = Flow()  # still considered via b1
            es.add(b1)
            # optimal

    try:
        model = Model(es)  # will fail for case 7
        solver_results = model.solve("cbc")
        termination = solver_results["Solver"][0]["Termination condition"]
    except:
        termination = "unable to build"

    print(variation, ": ", termination)


for variation in range(8):
    solve_case(variation)

[p-snft]

I tried to check indirectly by creating a networkx graph. However, that package implicitly adds nodes an edge connects with. Thus, I created oemof/oemof-network#51.