Why can't this simple oemof-model be optimized?

I also asked this question at stackexchange, but probably it’s in better hands here.

I would like to start with a very simple model. It consists only of a source (SRC, fixed), a sink (SNK, fixed) and a storage (STO, half filled). All three elements are connected to the same bus (BUS).
I would expect if SRC>SNK the STO is charged and if SRC<SNK the STO is discharged.
In both cases, nothing happens to STO and I get the following error message:

ERROR:root:Optimization failed with status ok and terminal condition infeasible

If SRC=SNK no errors are shown.

My code looks like:

import oemof.solph as solph
from oemof.outputlib import processing

es = solph.EnergySystem()    
BUS = solph.Bus(label='BUS')

SRC = solph.Source(label='SRC',
                      outputs={BUS: solph.Flow(actual_value=0.8,
                                                  nominal_value=100,
                                                  fixed=True)})
STO = solph.components.GenericStorage(label='STO',
                                         nominal_capacity=100000,
                                         initial_capacity=0.5,
                                         inputs={BUS: solph.Flow()},
                                         outputs={BUS: solph.Flow()})
SNK = solph.Sink(label='SNK',
                    inputs={BUS: solph.Flow(fixed=True,
                                               actual_value=0.5,
                                               nominal_value=100)})

es.add(BUS, SRC, STO, SNK)
om = solph.Model(es)
om.solve(solver='cbc', solve_kwargs={'tee': True})
print(processing.results(om))

What am I missing? :thinking:

obviously my fault!
the initial_capacity of the STO is also the final capacity. so the capacity can’t change if optimizing only one time step.

1 Like

You are right! the Generic_Storage capacity rule is one reason why single-timestep optimization doesn’t work at the moment. A workaround here is to loop over two time steps, neglect the results for the second time step, and update the initial_capacity. I have once implemented this in a Rolling Horizon analysis,

This approach is not really quick! If you need further information here, let me know!

As quite some users necessitate single-time-step optimization, the developer group will put it on our agenda for the next developer meeting in October.

1 Like

Does it work if you optimize multiple time-steps?

thx for your reply.
unfortunately in my simple model the suggested workaround won’t fix the problem.
the energy has to be balance! since the storages initial capacity is the final capacity, the energy provided by the fixed source must equal the demand of the fixed sink. the horizon or number of time steps won’t change anything i fear.
an other idea: is there a way to extract the objective function and the constraints of the lp inside the model?

Hi Niels,

as you noticed, the optimization for one time step does not work with components holding intertemporal constraints such as a storage.

I don’t know your exact problem but assuming that you are repeating the solving process for one time step, it might be a workarround to deal with intertemporal constraints in the repetition (iteration, …) manually.

For one timestep you could create a “storage bus” as follows which withdraws/releases energy from/to the electrical bus:

import oemof.solph as solph
from oemof.outputlib import processing, views

es = solph.EnergySystem()
BUS_EL = solph.Bus(label='BUS_EL')
BUS_STO = solph.Bus(label='BUS_STO')

SRC = solph.Source(
    label='SRC', outputs={BUS_EL: solph.Flow(
        actual_value=0.8, nominal_value=100, fixed=True)})

SNK = solph.Sink(
    label='SNK', inputs={BUS_EL: solph.Flow(
        fixed=True, actual_value=0.5, nominal_value=100)})

STO_IN = solph.Transformer(
    label='STO_IN',
    inputs={BUS_EL: solph.Flow(nominal_value=100)},
    outputs={BUS_STO: solph.Flow()})

STO_OUT = solph.Transformer(
    label='STO_OUT',
    inputs={BUS_STO: solph.Flow()},
    outputs={BUS_EL: solph.Flow(nominal_value=100)})

STO_SRC = solph.Source(label='STO_SRC', outputs={BUS_STO: solph.Flow()})
STO_SNK = solph.Source(label='STO_SNK', inputs={BUS_STO: solph.Flow()})

es.add(BUS_EL, BUS_STO, SRC, SNK, STO_IN, STO_OUT, STO_SRC, STO_SNK)

om = solph.Model(es)
om.solve(solver='cbc', solve_kwargs={'tee': True})

results = processing.results(om)

data = views.node(results, BUS_EL)

print(data['sequences'].head())

You can then deal with the intertemporal constrains manually e.g. by dealing with the storage level in a variable and creating a logic arround it. Or maybe even by creating a storage class which has a method that takes the optimization model as an argument in some method which is called in every iteration, …

Anyway, I am assuming that your real problem is quite big and optimization is the right method :wink:

Would this be a solution?

Cheers,
Cord

2 Likes

Of course I meant to model transformers with the bus. I have adapted my example accordingly…

thx!
this is a feasible solution to my problem, which is indeed more complex than the shown sniplet. :+1:

Thanks @CKaldemeyer for this solution. It is also very useful if you have a storage with a defined initial capacity (e.g. 0) and want to optimise the usage of the storage but not forcing the storage value back to the initial value in the last timestep. In this way I solved the problem for my long term storage hydrogen pressure tank.

I am glad you found a workaround.

Just for your interest. There is an open issue to make it possible to decouple the initial capacity and the capacity of the last time step in the future. It is just a lack of time that it is still open, sorry.

Certainly any help is welcome because I will not have time until December/January.