Understanding the OffsetTransformer

oemof

#1

Hello,

I am trying to understand oemof with one small unit commitment example, which I have implemented in pyomo a) as an abstract model and b) using blocks. There are only three time steps and three plants with respective economic dispatch and unit commitment constraints as well as some costs.

The following data is stored in a subfolder called input as csv files:
**********demand.csv ***********
T,Demand
1,150
2,500
3,400
************plant_params csv **********
J,Pmax,Pmin,RU,RD,FC,VC,SU,SD
1,370,50,200,300,200,4,800,20
2,200,80,100,150,280,5,720,12
3,140,40,100,100,240,6,200,40

The demand is self explaining. The parameters are in the column order: maximum/ minimum power output, ramp up/ down gradient per time step, fixed cost if in operation, variable cost, start up/ shut down cost.

At first I got an error, because I passed the OffsetTransformer coefficients as a dict, as it is presented in the class documentation. Then I realized that it should be a nested list according to the example included in the tests. Are these two different objects, which I don’t see so far, or is there maybe an inconsistency within the docs? Now I get the following:
WARNING:pyomo.core:Constant objective detected, replacing with a placeholder to prevent solver failure.
ERROR:root:Optimization failed with status ok and terminal condition infeasible

This is the code I am using for oemof:

import pandas as pd
from oemof.tools import helpers
import oemof.outputlib as outputlib
import oemof.solph as solph

demand = pd.read_csv('input/demand.csv')
plants = pd.read_csv('input/plant_params.csv', index_col=0)

solver = 'cbc'
energysystem = solph.EnergySystem(timeindex=demand.index)
bfuel = solph.Bus(label="fuel", variable_costs=1)
bel = solph.Bus(label="electricity")
energysystem.add(bfuel, bel)
energysystem.add(solph.Source(label='fuel_source',
                 outputs={bfuel: solph.Flow(nominal_value=10e3)}))

fixed_param = lambda val: [val for _ in demand.index]

for p in plants.index:
    plant = plants.loc[p]
    transf = solph.custom.OffsetTransformer(
        label="plant{}".format(plant.name),
        inputs={bfuel: solph.Flow(
                    nominal_value=plant['Pmax'],
                    min=fixed_param(plant['Pmax'] / plant['Pmin']), 
                    positive_gradient={'ub': plant['RU'], 'costs': 0},
                    negative_gradient={'ub': plant['RD'], 'costs': 0},
                    startup_costs=plant['SU'],
                    shutdown_costs=plant['SD'],
                    nonconvex=solph.NonConvex())},
        outputs={bel: solph.Flow()},
        coefficients=[fixed_param(plant['FC']),
                      fixed_param(plant['VC'])
                       ])
    energysystem.add(transf)


energysystem.add(solph.Sink(label='demand', inputs={bel: solph.Flow(
    actual_value=demand['Demand'], fixed=True, nominal_value=1)}))

model = solph.Model(energysystem)
model.solve(solver=solver)

I would be very happy, if anyone has an idea, what I do wrong.
Best Philipp


#2

Hi Philipp,

your example does not work for me with v.0.2.1 even if I fix the missing datetimeindex for the demand-dataframe. Which oemof-version do you have installed?

There’s also a test on the offsettransformer which should work with v.0.2.2 at which you might have a look.

Cheers,
Cord


#3

Hello Philipp,

the first rule for debugging is to minimise the problem to reduce sources of errors.

  1. Use only one transformer that can satisfy the demand
  2. Write the parameters directly into the class definition: startup_costs=5 instead of startup_costs=plant['SU']
  3. Try the component with the essential parameters first and add the others (gradients, startup_costs…) step by step.
  4. Write out the lp-file to check the problem:
    model = solph.Model(energysystem)
    model.write('/your/path/my.lp', io_options={'symbolic_solver_labels': True})
    model.solve(solver=solver)

Furthermore, you can avoid the warnings:

  1. Create a datetime index:
- energysystem = solph.EnergySystem(timeindex=demand.index)
+ date_time_index = pd.date_range('1/1/2012', periods=len(demand),freq='H')
+ energysystem = solph.EnergySystem(timeindex=date_time_index)
  1. Add the variable_costs to the Flow and not to the Bus
- bfuel = solph.Bus(label="fuel", variable_costs=1)
- energysystem.add(solph.Source(label='fuel_source',
-                  outputs={bfuel: solph.Flow(nominal_value=10e3)}))
+ bfuel = solph.Bus(label="fuel")
+ energysystem.add(solph.Source(label='fuel_source',
+                  outputs={bfuel: solph.Flow(nominal_value=10e3, variable_costs=1)}))

The component you use is still marked as experimental so it is poorly documented and may contain little bugs. But it should work in general. One little bug is that you have to pass lists a coefficients. Normally you should be able to pass a list or a scalar. Scalars are useful if all values are the same: just [5, 0.4] instead of [[5, 5, 5], [0.4, 0.4, 0.4, ]]. I will fix that but until this is published you have to use the second variant.

If it still does not work please post your new (small) example and the error message again.


#4

I was already about to answer but Uwe adressed all of the points :wink:

Here’s my (obsolete) answer:

I now had a look with the latest dev Branch. The documentation is wrong somehow and you are right that it is a nested list.

This error must be context-wise and related to your problem/setup.

You might want to check the LP-file using:

model.write('example.lp', io_options={'symbolic_solver_labels': True})

Best,
Cord


#5

Your min parameter is greater than one (!). The max value is one by default. The variable is limited by its max and min. The variable has to be greater than
min * nominal value and (!) smaller than max * nominal_value.

That’s why it is infeasible.

All these parameters (min, max,…) are relative to the nominal value. The code below works for me!

In your case the nominal value is something like the installed capacity. This value should not be a list. It is better to use min/max as lists to reduce the installed capacity.

Use shortage and excess variables to avoid infeasible problems (see Solver time seems to be to high for a small energy model using oemof (solph) with CBC solver ). It is easier to find errors if some results exist.

import pandas as pd
import oemof.solph as solph

demand = pd.read_csv('data/demand.csv')
plants = pd.read_csv('data/plant_params.csv', index_col=0)

date_time_index = pd.date_range('1/1/2012', periods=len(demand),
                                freq='H')

solver = 'cbc'
energysystem = solph.EnergySystem(timeindex=date_time_index)

bfuel = solph.Bus(label="fuel")
bel = solph.Bus(label="electricity")
energysystem.add(bfuel, bel)
energysystem.add(solph.Source(label='fuel_source',
                 outputs={bfuel: solph.Flow(nominal_value=10e3,
                                            variable_costs=1)}))

fixed_param = lambda val: [val for _ in demand.index]


for p in plants.index:
    plant = plants.loc[p]
    transf = solph.custom.OffsetTransformer(
        label="plant{}".format(plant.name),
        inputs={bfuel: solph.Flow(
                    nominal_value=400,
                    min=fixed_param(plant['Pmax'] / (plant['Pmin'] * 10)),
                    positive_gradient={'ub': plant['RU'], 'costs': 0},
                    negative_gradient={'ub': plant['RD'], 'costs': 0},
                    startup_costs=plant['SU'],
                    shutdown_costs=plant['SD'],
                    nonconvex=solph.NonConvex())},
        outputs={bel: solph.Flow()},
        coefficients=[[0, 0, 0], [1, 1, 1]])
    energysystem.add(transf)


energysystem.add(solph.Sink(label='demand', inputs={bel: solph.Flow(
    actual_value=demand['Demand'], fixed=True, nominal_value=1)}))

model = solph.Model(energysystem)
model.solve(solver=solver)

#6

Thanks a lot for the fast reply of you two!

@CKaldemeyer
I’m using v 0.2.2. I also played with the set of the periods. Internally, it is quickly converted to natural numbers anyway, I think. Uwe addressed the warning. Furthermore, I used the example from the test mentioned as I mentioned above.

@uwe.krien
Thank you for the valuable information! I haven’t checked yet, but it makes sense. Actually, it also felt a bit strange to use a slope > 1. Therefore, I tricked with the fuel price.
I already said to Cord on the Phone, that I intend on using oemof more often. So I would be glad to contribute somehow. In this case I could extend the docs with the min and max restrictions. I haven’t found these in OffsetTransformer nor in Transformer which it inherits from.


#7

You are very welcome :smiley:

The min/max attribute belongs to the Flow, therefore you have to look in the docs of the Flow class. Once you have understood the Flow attributes you can use them with all components (Sink, Source, Transformer,…).

Please check your example and set this post to “solved” if everything works fine.


#8

@uwe.krien
I definitely should fully understand the Flow class as well as its interpretation. It seems to be pivotal.
My example is now running, but the result is not as expected. This is probably due to the fact, that the intended fixed and variable costs are not considered anymore in the example (coefficients of 0 and 1 instead of ‘FC’ and ‘VC’). If I get it running properly, I will post the solution. After that, I close the post.


#9

I tricked around and got it working. It actually looks a bit nasty, because the example is taken from an excercise and contains a fixed and variable cost approach, which does not (to my current knowledge) fit to the reality or the oemof philosophy. Therefore, some mangling to convert the cost function to a transformer input-output is necessary. Here is the solution:

import pandas as pd
from oemof.tools import helpers
import oemof.outputlib as outputlib
import oemof.solph as solph

demand = pd.read_csv('input/demand.csv')
plants = pd.read_csv('input/plant_params.csv', index_col=0)

fixed_param = lambda val: [val for _ in demand.index]
linear_cost_inverse = lambda fc, vc: (-fc / vc, 1 / vc)

def plant_coefficients(s_plant, n_periods=len(demand)):
    """ return OffsetTransformer coefficients from plant """
    fc, vc = s_plant['FC'], s_plant['VC']
    n, m = linear_cost_inverse(fc, vc)
    return [[n for _ in range(n_periods)], [m for _ in range(n_periods)]]


min_fuel = lambda s_plant: s_plant['FC'] + s_plant['Pmin'] * s_plant['VC']
max_fuel = lambda s_plant: s_plant['FC'] + s_plant['Pmax'] * s_plant['VC']

assert linear_cost_inverse(200, 4) == (-50, 0.25)
assert plant_coefficients(plants.loc[1]) == [[-50, -50, -50],
                                            [0.25, 0.25, 0.25]]
assert min_fuel(plants.loc[1]) == 400
assert max_fuel(plants.loc[1]) == 1680
solver = 'cbc'
date_time_index = pd.date_range('2012-01-01', periods=len(demand.index),
                                freq='H')

energysystem = solph.EnergySystem(timeindex=date_time_index)
bfuel = solph.Bus(label="fuel")
bel = solph.Bus(label="electricity")
energysystem.add(bfuel, bel)
energysystem.add(solph.Source(label='source',
                 outputs={bfuel: solph.Flow(nominal_value=10e3,
                                            variable_costs=1)}))

for p in plants.index:
    plant = plants.loc[p]
    transf = solph.custom.OffsetTransformer(
        label="plant{}".format(plant.name),
        inputs={bfuel: solph.Flow(
                    nominal_value=max_fuel(plant),
                    min=fixed_param(min_fuel(plant)/max_fuel(plant)),
                    startup_costs=plant['SU'],
                    shutdown_costs=plant['SD'],
                    nonconvex=solph.NonConvex())},
        outputs={bel: solph.Flow(nominal_value=plant['Pmax'],
                    positive_gradient={'ub': plant['RU'], 'costs': 0},
                    negative_gradient={'ub': plant['RD'], 'costs': 0},
                                )},
        coefficients=plant_coefficients(plant))
    energysystem.add(transf)


energysystem.add(solph.Sink(label='demand', inputs={bel: solph.Flow(
    actual_value=demand['Demand'], fixed=True, nominal_value=1)}))

model = solph.Model(energysystem)
model.solve(solver=solver)
energysystem.results['main'] = outputlib.processing.results(model)
energysystem.results['meta'] = outputlib.processing.meta_results(model)
results = energysystem.results['main']

assert model.objective.expr() == 7572
# current output [[150, 370, 360], [0, 130, 0], [0, 0, 40]]
assert expected_plant_dispatch == [[150, 350, 360], [0, 100, 0], [0, 50, 40]]

Some questions remain:

  1. Is there an error with the ramp up gradient? The included data restricts it to 200 for plant one. The feasible output exhibits a gradient of 220 from time step one to two.
  2. Is there a reason that the NonConvex Flow is implemented at the input flow? It seems more natural or it would be nice to define these plant properties to the output.
  3. Plants might have fixed cost per hour. The current implementation puts all the emphesis on the input/fuel use. It somehow could be pressed into this form. But it actually is very cumbersome with varying fuel prices. I think the easiest ist to create an additional term in the objective function. Is this already under development? Otherwise, I would put work load into this, because it is crucial for many practical applications.

Cheers


#10
  1. All these values are relative to the nominal value (see my post above), so please check if:
    nominal_value * positive_gradient = 200.

  2. Normally it does not matter if you define the values for the input or the output flow but in the OffsetTransformer the non-convex flow has to be in the input flow. I don’t know why.
    Maybe @CKaldemeyer does know more about it.

  3. This is already under development in PR509. Ask @sar_b for more information.


#11

I don’t have time to reproduce your example but will answer on your second question:

Long answer:

I guess the term “intuitive” is very subjectively perceived. A while ago, I have implemented the component quite quickly because I needed it in this way in a project to model heat pumps. At the moment, the component is under custom (which could also be named “quite experimental”) and the selection of the input flow for the binary was initially due to the idea of extending the component to have multiple output flows with single conversion factors. Then, “muting” the output flows would be easy this way.

But of course, the component could be described vice versa or the binary flow could be detected and the equations selected depending on the equation side. Or we could try no enable a n:m relation… Feel free to add this functionality using a PR! :wink:

Short answer:

No, there is not particular reason and I needed it that way when I added it to the code base.

Best wishes,
Cord


#12

If you need this “offset behaviour” for back-pressure turbines, you could also use the GenericCHP component.

An example can be found here and a describing article here. Note that I have adapted the approach for the flue gas losses (which are passed as a param) as described in the docstring.

Best,
Cord


#13

Thank you @uwe.krien for the explaination. With the gradient relative to the nominal_value, the dispatch now seems to fit. I will have to check for the objective later. Anyway, I think I can close it now.

@CKaldemeyer: Thank you for the historical story to it. Then, things often make more sense. As said, I would love to contribute. Unfortunately, I do have some things to finish, before being able to return to oemof. :frowning: I do fully agree that intuitive is very subjective. In this case I connected it to the “common” wording, which technicians tend to. With power plants people talk more often about minimal power output than minimal fuel use, when characterizing a plant, I would say.

I had already put an eye on your CHP formulation :wink:
Best whishes,
Philipp