Initial investment costs

[henhuy]

I'm trying to implement initial costs for a PV component which holds an investment optimization. For example, add 100€ to PV costs if (and only if) PV is built: 0 kW -> 0€ x kW -> 100€ + x kW * 2 €/kW

As far as I can see, there is not yet an option for that in oemof. I tried to implement an additional InitialInvestmentFlow but as this leads to RuntimeError: Cannot write legal LP file. Objective 'objective' has nonlinear terms that are not quadratic. Is it even possible to implement it? Some hints would be nice!

See my code example below. Problematic line is:

initial_costs += (
    (sum(m.flow[i, o, t] for t in m.TIMESTEPS) > 0) *
    m.flows[i, o].investment.initial_costs
)

If I remove "greater than" part, some costs are added (not what I want, but only to check if custom flow works)

Example script:

"""
Test scratch to implement pv initial invest
"""

import pandas
from oemof.solph import EnergySystem, Bus, Flow, Sink, Source, Investment
from oemof.solph import Model
from oemof import outputlib

from pyomo.core import Set, Expression
from pyomo.core.base.block import SimpleBlock
import oemof.groupings as groupings

class InitialInvestmentFlow(SimpleBlock):
    def _create(self, group=None):
        if group is None:
            return None

        self.FLOWS = Set(initialize=[(g[0], g[1]) for g in group])

    def _objective_expression(self):
        if not hasattr(self, 'FLOWS'):
            return 0

        m = self.parent_block()
        initial_costs = 0

        for i, o in self.FLOWS:
            if m.flows[i, o].investment.initial_costs is not None:
                initial_costs += (
                    (sum(m.flow[i, o, t] for t in m.TIMESTEPS) > 0) *
                    m.flows[i, o].investment.initial_costs
                )
            else:
                raise ValueError("Missing value for fixed investment costs!")

        self.initial_costs = Expression(expr=initial_costs)
        return initial_costs

def _initial_investment_grouping(stf):
    if hasattr(stf[2], 'investment'):
        if stf[2].investment is not None:
            return True
    else:
        return False

initial_investment_flow_grouping = groupings.FlowsWithNodes(
    constant_key=InitialInvestmentFlow,
    filter=_initial_investment_grouping
)

energysystem = EnergySystem(
    timeindex=pandas.date_range('2019-01-01', periods=2, freq='H'),
    groupings=[initial_investment_flow_grouping]
)
bus = Bus()
net = Source(
    label="Net",
    outputs={
        bus: Flow(nominal_value=90)
    }
)
investment = Investment(ep_costs=2)
investment.initial_costs = 100
pv = Source(
    label="PV",
    outputs={
        bus: Flow(
            actual_value=[0.2, 0.2],
            fixed=True,
            investment=investment,
        )
    }
)
demand = Sink(
    label="demand",
    inputs={
        bus: Flow(
            nominal_value=1,
            actual_value=[100, 100],
            fixed=True
        )
    }
)
excess = Sink(
    label="excess",
    inputs={bus: Flow()}
)

energysystem.add(bus, net, pv, demand, excess)
om = Model(energysystem=energysystem, constraint_groups=[InitialInvestmentFlow])
om.solve(
    solver="cbc",
    solve_kwargs={'tee': True, 'keepfiles': True}
)

results = outputlib.processing.results(om)
# import pdb;pdb.set_trace()
print(results)

[simnh]

Hey, as far as I see the problem, the only way would be to use binary variables. Adding

x <= y * investment_max 
initial_cost = y * 100

with y being binary {0,1}, x the investment

So y must be 1 if x is supposed to be anything then zero and otherwise 0. Adding initial cost only in the case of investment.

investment_max can also be an arbitrary high number, but close to the highest value you would expect for x to avoid numerical instability of the problem due to differences the matrix (Big-M-constraint).

Actually you could take a look at the offset transformer, which is mathematically a similar if not same problem (if I remember correctly). I was on holiday quite a while an my head is not yet working a 100%, but I think it should be right and unfortunately also the only solution to your problem (increasing run time significantly). Correct me if you find an error, but you will need a binary variable for investment yes/no for sure.

[henhuy]

Hey @simnh , thank you for your hint! I will look into that. Nice

[joroeder]

hey, some time ago, we were discussing this topic within our project as well. So, thanks for opening this issue!

We tried and implemented some stuff, but did not manage to propose it here on github so far. And we realized the same problem, @simnh was mentioning:

investment_max can also be an arbitrary high number, but close to the highest value you would expect for x to avoid numerical instability of the problem due to differences the matrix (Big-M-constraint).

I will open a PR for proposing and discussing our implementation. I think that's a nice feature, which should not be missed in oemof :)