Modeling PV-Feed-In



Hey everyone,

I am trying to model a system which supplies its electrical energy demand either out of grid energy which causes costs for the operator or supplies the demand out of renewable sources e.g. pv. Excess pv-energy is fed into the grid which causes negative costs or income for the operator. This is modeled by a two-way-transformer see sketch below:

The simulations result is an energy flow in both directions at the same timestep which I want to avoid. So my attempt was to formulate a constraint using the BinaryFlow-Class:

my_block = environ.Block()

def connect_flow_rule(m,t):

       expr = (om.BinaryFlow.status[load, Bus1, t] *  om.flow[Bus1, feed,t] == 0

return expr

my_block.flow_connect_constr = environ.Constraint(om.TIMESTEPS, rule = connect_flow_rule)
om.add_component('ConnectFlow', my_block)

This causes a nonlinear problem similar to the topic in : Add constraints to Binary Flow
I guess there must be another formulation of the constraint but I can’t see it. Maybe someone has already had the same problem and can give advice!

Thanks in advance!


Hello Bastian. Didn’t this question come up at the Frankfurt workshop? I have three responses. First, does the solver actually choose the loop flows you suggest, that is, are you 100% certain there is a problem? I suspect the solver does not do this but I cannot cite a reliable source for that view. But neither have I seen a minimum cost flow problem algorithm that adds the cost-less loop flows you indicate simply because “it can”. Second, you could add minuscule costs or losses to the flows in question: sufficient to influence the solver but not so large as to make a material difference to the results. Or you can go down the binary route and limit the flow to one direction only. This is computationally expensive though. (As an aside, I once implemented a nodal pricing algorithm using this third option. Not only is it slow, it also yields faulty marginal costs. As Daniel Huppmann explains here. I duly filed a bug report to myself!) HTH, Robbie. PS: nice diagram, makes the question easier to visualize.


Hello Robbie,
thank you very much for your reply and sorry for my late response. I havent’t been to the workshop in Frankfurt…
Anyways, yes I am definitely sure that the solver choose the “loop energy flow”. I could see it at the simulations results sheet which I am extracting as a csv over all timesteps.
The minuscule costs you said reduced the loop energy flow but couldn’t totally avoid it. But that got me on the track… I modeled “the grid” with a unspecified source/sink like for example commodity sources. I got rid of the grids sink/source and changed the bus attribute to “balanced=FALSE” which disables the conservation of energy at the bus balance. The “loop energy flows” disappeared subsequently. Maybe a thing that could go to oemof-documentation…
Anyway that problem is solved, thanks a lot and also for the paper…
Best regards, Bastian


Hello Bastian. That’s good you have a solution. But intuitive results do not necessarily imply a correct model! To be blunt, I am more inclined to trust the solver than your coding? Just my reaction from a considerable distance. Robbie :neutral_face:


I agree with you but I don’t want to trust the solvers results due to a “loop energy flow” because it is physically not possible. In my eyes it is more just a theoretical result which is probably more optimized but cannot happen in reality. I checked the solutions results extensively and it should reproduce the reality in a physical correct way. Bastian


Hi Bastian,

the product of a continous and integer variable can be resolved using the Big M method.

Your case (if y>0: x=0; if x>0: y=0) can also be resolved using this method which results in a classic MILP problem with 4 additional constraints. There are several techniques to model boolean logic, etc. which can be found in the AIMMS modelling guide or (if you are familiar which the German language, which I presume) in the following book:

Kallrath, Gemischt-ganzzahlige Optimierung: Modellierung in der Praxis (2. Auflage), Springer Vieweg, 2012

I could also send you my master thesis where I have modelled exactly your problem case applied to another topic.

Good luck!


Unfortunately, there is no LP formulation for this which I have already asked on stackoverflow because I was hoping for some modelling guru tweaking constraints :wink:




Hello together!
One of my master students has worked on a very similar problem and we formulated it as a MILP.
Kind regards, Martin


I guess this should do it:

feed(t) <= feed_max * y(t)
load(t) <= load_max * (1-y(t))

where feed/load(t) are respective continous decision variables (oemof flows), feed/load_max are parameters/ upper bounds (nominal capacities of oemof flows) and y(t) is a binary decision variable.

Good luck!


I guess this topic could be marked as resolved. Bastian, will you do it if you have tested our proposed solutions?


Hello @CKaldemeyer and @martin.klein. Are you sure this solution meets the requirements of the original question? As I understand it, the load has preferential access to the PV facility because, I guess, they own the installation and because it has a SRMC (short-run marginal cost) of supply of zero. But didn’t the owner also want to sell surplus PV generation back to the grid operator. In which case, toggling the PV facility between load (supplying the owner) mode and feed-in (back to the grid) mode does not deal with the excess output case. Or am I missing something?

My other comment is that the use of constraint programming in this way is only necessary if the research question involves PV owners and their actions and incentives. If not, then why not just treat every load and generator element as part of the one undifferentiated problem? While acknowledging that least cost for the system is not necessarily least cost for the PV owner. Or alternatively, model the incentives of every actor and not just one subset, in this case PV owners. Robbie


Hello together,

@robbie.morrison, the research question in that case is in fact based on the PV owners point of view.

Thank you @CKaldemeyer for your input. I solved the system with a constraint as followed:

my_block_1 = environ.Block()

def storage_lock_IN_rule(m,t):
      expr = om.BinaryFlow.status[b_el, el_storage, t] * el_storage.nominal_capacity = om.flow[b_el, el_storage,t]
      return expr

my_block_1.storage_lock_IN_constr = environ.Constraint(om.TIMESTEPS, rule = storage_lock_IN_rule)
om.add_component('StorageLockIN', my_block_1)

my_block_2 = environ.Block()

def storage_lock_OUT_rule(m,t):
      expr = (1 - om.BinaryFlow.status[b_el, el_storage, t]) * el_storage.nominal_capacity = om.flow[el_storage, b_el, t]
      return expr

my_block_2.storage_lock_OUT_constr = environ.Constraint(om.TIMESTEPS, rule = storage_lock_OUT_rule)
om.add_component('StorageLockOUT', my_block_2)

This example is about avoiding loop energy flow in storage’s. But the approach of @CKaldemeyer could be used generally to avoid loop energy flows if the research question requires constraint programming. But indeed the use of BinaryFlows results in a longer simulation time.


Hello all. Let me be philosophical for a moment. Under the programming doctrine of separation of concerns and noting that engineering performance and (sub)system dispatch are fully orthogonal concepts, I prefer that system control decisions not be embedded in plant definitions. I have implemented this design approach so it can work and it will reduce both the class count (and maintenance overhead) and improve the design space significantly.

In passing, storage entered this discussion quite late on: use of storage and other cross-temporal issues like DSM add to the problem complexity, but more importantly to the software design challenge. Robbie