Constraint setting flow value to Zero

Q&A
1

Hello openmod world,
as I’m new to the community don’t hesitate to give me feedback on my way of asking questions in this forum.
I’m trying to model the variable revenue for electricity generation of a block heat and power plant (BHKW) connected with the investment method. Due to variable fundings for electricity generation of a BHKW depending on full load hours (3500h), it is necessary to vary the variable output costs of the component after reaching this amount. An idea so far is to define two electricity busses to the component with different variable output costs an adding a constraint, or two.
First to calculate the full load hours after every timestep based on investment capacity.
Second to set the flow value of the low revenue bus to zero until reaching the amount of funded full load hours and afterwards setting the high revenue bus flow to zero.

Does anybody has an idea how to implement this behavior of a transformator?
Or is there already an example how to do it?
Is it even possible to couple the investment optimization to two output buses for this case?

I’m looking forward to any suggestions.
Kind regards and have a nice day!

2

Hi Mathias,

I set up an example that uses two auxiliary transformers, an auxiliary electricity bus and an addional constraint to model the desired behavior. The electricity from the CHP can go through both of the transformers to get to the demand sink. Both have a conversion factor of 1. One provides no revenue (variable_costs=0) while the other transformer provides an addional income (ger.: “KWK-Vergütung”) for a certain amount of hours. The electricity price itself (base tarif, spot or whatever you are considering) is still set in the demand sink. The limitation then is given via an addional constraint that sets the summed flow through the “high revenue”-transformer equal to the amount of electricity produced during the full load hours (nominal value of the CHP*full load hours).

Best,
Jakob

import pandas as pd
from oemof import solph
import oemof.outputlib as outputlib
import pyomo.environ as po

date_time_index = pd.date_range('1/1/2012', periods=50, freq='H')
energysystem = solph.EnergySystem(timeindex=date_time_index)
b_gas = solph.Bus(label='gas')
b_el = solph.Bus(label='electricity')
# Auxiliary electricity bus to realize the additional constraint
b_aux_el = solph.Bus(label='aux_el')
b_heat = solph.Bus(label='heat')
energysystem.add(b_gas, b_el, b_heat, b_aux_el)

energysystem.add(solph.Source(label='gassupply',
                              outputs={b_gas: solph.Flow(variable_costs=1)}))

energysystem.add(solph.Sink(label='grid',
                            inputs={b_el: solph.Flow(variable_costs=-2)}))

energysystem.add(solph.Sink(label='DH',
                            inputs={b_heat: solph.Flow(variable_costs=-1,
                                                       nominal_value=5,
                                                       actual_value=1,
                                                       fixed=True)}))

energysystem.add(solph.Transformer(
    label='CHP',
    inputs={b_gas: solph.Flow()},
    outputs={
        b_heat: solph.Flow(),
        b_aux_el: solph.Flow(
            investment=solph.Investment(ep_costs=5))}, 
    conversion_factors={b_aux_el: 0.3, b_heat: 0.4}))

energysystem.add(solph.Transformer(
    label='aux_trafo_01',
    inputs={b_aux_el: solph.Flow()},
    outputs={b_el: solph.Flow()},
    conversion_factor=1,
    variable_costs=0))

# Full load hours (flh): determines the amount of electricity with high_revenue
flh = 10
# KWK-Vergütung (additional revenue for electricity from CHP)
addition_chp_revenue = -1
# Auxiliary component that provides additional income for electricity from
# the CHP for a certain amount of electricity (german:"KWK-Vergütung").
# The 'summed_max' attribute of this component will be used in the
# additional constraint.
energysystem.add(solph.Transformer(
    label='aux_trafo_02',
    inputs={b_aux_el: solph.Flow()},
    outputs={b_el: solph.Flow()},
    conversion_factor=1,
    variable_costs=addition_chp_revenue))

om = solph.Model(energysystem)
chp = energysystem.groups['CHP']
high_revenue_trafo = energysystem.groups['aux_trafo_02']
# grid = energysystem.groups['grid']

# Add an additional constraint that limits the amount of electricity that can
# go through aux_trafo_02
myconstrains = po.Block()
om.add_component('MyBlock', myconstrains)
myconstrains.solar_constr = po.Constraint(
    expr=((sum(om.flow[high_revenue_trafo, b_el, t] for t in om.TIMESTEPS))
          <= flh*om.InvestmentFlow.invest[chp, b_aux_el]))

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

energysystem.results['main'] = outputlib.processing.results(om)
energysystem.results['meta'] = outputlib.processing.meta_results(om)
results = energysystem.results['main']

nominal_invest_CHP = (outputlib.views.node(results, 'CHP')[
    'scalars'][(('CHP', 'aux_el'), 'invest')])
electricity_bus = outputlib.views.node(results, 'electricity')
heat_bus = outputlib.views.node(results, 'heat')
print('********* Main results *********')
print(electricity_bus['sequences'].sum(axis=0))
print(heat_bus['sequences'].sum(axis=0))
print("")
print("Capacity (electr.) of CHP (nominal_value):", nominal_invest_CHP)
print("Total amount of full load hours:",
      electricity_bus['sequences'].values.sum()/nominal_invest_CHP, "h")
print("Amount of full load hours with higher revenue: ", flh, "h")
3

HiJakob,

thanks a lot for your quick response! Testing your system I realized, that all electricity generated by the chp is only distributed by the aux_trafo_01 without additional revenue (3.75 kW, 187.5 kWh). The flow from aux_trafo_02 to grid is zero. I tried to put the variable_costs definition into the solph.Flow()-argument inside the output definition of these transformers. This resulted in a higher nominal value of the chp (18.75 kW) because of the invest optimization and all generated electricity was distributed by aux_trafo_02 (187.5 kWh). Setting a maximum constraint for the investment solved this problem and the electricity was distributed among the two aux_trafos (3.75 kW, 150 kWh, 37.5 kWh). But this is not the correct workaround for the optimization i think. Furthermore the calculation of the total full load hours can not be calculated with .values.sum() of the electricity_bus cause it sums the values of input and output of the bus I think. This results in full load hours value twice as high as correct. I’ll see if I find a way to adress the variable revenue to the correct flow, that the invest optimization works properly.
If anybody is interested in this topic, or has additional suggestions, don’t hesitate to participate :slight_smile:
Kind regards,
Mathias

4

Hi Jakob,
with some help of Uwe I could solve the problem of the different component investment sizes. Some additional constraints connecting the invest variables of the trafos were necessary. I’ll post the adapted code below. I hope this works…
Best regards!

import pandas as pd
from oemof import solph
import oemof.outputlib as outputlib
import pyomo.environ as po

date_time_index = pd.date_range('1/1/2012', periods=50, freq='H')
energysystem = solph.EnergySystem(timeindex=date_time_index)
b_gas = solph.Bus(label='gas')
b_el = solph.Bus(label='electricity')
# Auxiliary electricity bus to realize the additional constraint
b_aux_el = solph.Bus(label='aux_el')
b_heat = solph.Bus(label='heat')
energysystem.add(b_gas, b_el, b_heat, b_aux_el)

energysystem.add(solph.Source(label='gassupply',
                              outputs={b_gas: solph.Flow(variable_costs=1)}))

energysystem.add(solph.Sink(label='grid',
                            inputs={b_el: solph.Flow(variable_costs=-2)}))

energysystem.add(solph.Sink(label='DH',
                            inputs={b_heat: solph.Flow(variable_costs=-1,
                                                       nominal_value=5,
                                                       actual_value=1,
                                                       fixed=True)}))

energysystem.add(solph.Transformer(
    label='CHP',
    inputs={b_gas: solph.Flow()},
    outputs={
        b_heat: solph.Flow(),
        b_aux_el: solph.Flow(variable_costs=-1,
                             investment=solph.Investment(ep_costs=10))},# maximum=3.75))},
    conversion_factors={b_aux_el: 0.3, b_heat: 0.4}))

# Full load hours (flh): determines the amount of electricity with high_revenue
flh = 10
# KWK-Vergütung (additional revenue for electricity from CHP)
addition_chp_revenue = -5

energysystem.add(solph.Transformer(
    label='aux_trafo_01',
    inputs={b_aux_el: solph.Flow()},
    outputs={b_el: solph.Flow(investment=solph.Investment(ep_costs=0))},
    variable_costs=0,
    conversion_factor=1
))

# Auxiliary component that provides additional income for electricity from
# the CHP for a certain amount of electricity (german:"KWK-Vergütung").
# The 'summed_max' attribute of this component will be used in the
# additional constraint.
energysystem.add(solph.Transformer(
    label='aux_trafo_02',
    inputs={b_aux_el: solph.Flow()},
    outputs={b_el: solph.Flow(investment=solph.Investment(ep_costs=0))},
    variable_costs=addition_chp_revenue,
    conversion_factor=1
))

om = solph.Model(energysystem)
chp = energysystem.groups['CHP']
low_revenue_trafo = energysystem.groups['aux_trafo_01']
high_revenue_trafo = energysystem.groups['aux_trafo_02']
# grid = energysystem.groups['grid']

# Add an additional constraint that limits the amount of electricity that can
# go through aux_trafo_02
myconstrains = po.Block()
om.add_component('MyBlock', myconstrains)
myconstrains.solar_constr = po.Constraint(
    expr=((sum(om.flow[high_revenue_trafo, b_el, t] for t in om.TIMESTEPS))
          <= flh*om.InvestmentFlow.invest[chp, b_aux_el]))

# Add additional constraints to connect investment variables
# of CHP, aux_trafo_01 and aux_trafo_02
solph.constraints.equate_variables(
    om,
    om.InvestmentFlow.invest[chp, b_aux_el],
    om.InvestmentFlow.invest[high_revenue_trafo, b_el])
solph.constraints.equate_variables(
    om,
    om.InvestmentFlow.invest[low_revenue_trafo, b_el],
    om.InvestmentFlow.invest[high_revenue_trafo, b_el]),
solph.constraints.equate_variables(
    om,
    om.InvestmentFlow.invest[chp, b_aux_el],
    om.InvestmentFlow.invest[low_revenue_trafo, b_el])

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

energysystem.results['main'] = outputlib.processing.results(om)
energysystem.results['meta'] = outputlib.processing.meta_results(om)
results = energysystem.results['main']

nominal_invest_CHP = (outputlib.views.node(results, 'CHP')[
    'scalars'][(('CHP', 'aux_el'), 'invest')])
electricity_bus = outputlib.views.node(results, 'electricity')
aux_bus = outputlib.views.node(results, 'aux_el')
heat_bus = outputlib.views.node(results, 'heat')
CHP = outputlib.views.node(results, 'CHP')
print('********* Main results *********')
print('seq_sum electricity_bus\n', electricity_bus['sequences'].sum(axis=0))
print('seq_max electricity_bus\n', electricity_bus['sequences'].max())
print('seq_sum aux_bus\n', aux_bus['sequences'].sum(axis=0))
print('seq_max aux_bus\n', aux_bus['sequences'].max())
print('seq_sum heat_bus\n', heat_bus['sequences'].sum(axis=0))
print("")
print("Capacity (electr.) of CHP (nominal_value):", nominal_invest_CHP)
print("Total amount of full load hours:",
      aux_bus['sequences'].sum(axis=0)/nominal_invest_CHP, "h")
print("Amount of full load hours with higher revenue: ", flh, "h")
1 Like
5

There is a little mistake in the definition of the aux_trafo_02-transformer:
variable_costs are a attribute for a flow, not for a transfomer. So in the example of the code they are ignored, which leads to wrong results.

1 Like
6

@PyMoF Please mark your solution as “solution”.

You can also edit your solution according to @Frapl.

Thank you both :slight_smile: