Hello everyone,
@EkaterinaZol and I are trying to do a stochastic optimization by sampling scenarios. Essentially we need several copies of the energy system with some components varying randomly, and other components identified across scenarios. At the last OEMOF developer meeting the strategy suggested was to use equate_variables, unfortunately that is not straightforward and we also wanted to look at another potential strategy.
Let me set up a toy problem (jupyter notebook is attached).
First we take some randomized demand, and some synthetic solar data:
number_of_time_steps = 24 * 3
date_time_index = pd.date_range("1/1/2012", periods=number_of_time_steps, freq="H")
solar_data = [np.sin(x.hour/24. * 2. * np.pi) if np.sin(x.hour/24. * 2. * np.pi) > 0. else 0. for x in date_time_index]
demand_data = [np.sin(x.hour/24. * 2. * np.pi) * 0.3 + 0.5 for x in date_time_index] + 0.05 * np.random.randn(len(date_time_index))
demand_data_2 = [np.sin(x.hour/24. * 2. * np.pi) * 0.3 + 0.5 for x in date_time_index] + 0.05 * np.random.randn(len(date_time_index))
Next we want to define the energy system that we want to sample scenarios with:
def scenario(tag, solar_data, dem_data):
bel = solph.Bus(label="electricity" + tag)
gsi = solph.Sink(
label="grid_sink"+ tag,
inputs={bel: solph.Flow(nominal_value=5000, variable_costs=0.1)},
)
gso = solph.Source(
label="grid_source"+ tag,
outputs={bel: solph.Flow(nominal_value=5000, variable_costs=4)},
)
dem = solph.Sink(
label="demand"+ tag,
inputs={bel: solph.Flow(fix=dem_data, nominal_value=10)},
)
pvi = solph.Source(
label="pv"+ tag,
outputs={bel: solph.Flow(fix=solar_data, investment=solph.Investment(ep_cost=0.1))},
)
return (bel, gsi, gso, dem, pvi)
Now we can make two copies of the scenario with the randomized demand:
es_2 = solph.EnergySystem(timeindex=date_time_index)
es_2.add(*scenario("_1", solar_data, demand_data))
es_2.add(*scenario("_2", solar_data, demand_data_2))
This creates two entirely independent scenarios. What we want, though is for the Investment and the grid_source flow (which we think of as purchased days ahead) to be the same across scenarios.
For the investment we ca do this with equate_variables:
m2 = solph.Model(es_2)
solph.constraints.equate_variables(m2,
m2.InvestmentFlow.invest[es_2.groups['pv_1'], es_2.groups['electricity_1']],
m2.InvestmentFlow.invest[es_2.groups['pv_2'], es_2.groups['electricity_2']])
This is copied pretty much from the documentation. The problem is that we don’t understand the underlying logic of oemof well enough to actually understand where we can find the variables we need to identify generically. Concretely how would we equate flow variables like the grid_source output? Further, it seems that this way we have two variables and we then introduce an equality among them. While fine for Investment variables if we have a time line this would lead to a lot of superfluous variables.
A logical alternative
What we would really like to do, but have no idea if it’s feasible, would be something along these lines:
# Doesn't work!
def scenario_broken(tag, solar_data, dem_data, pv_invest, source_flow):
bel = solph.Bus(label="electricity" + tag)
gsi = solph.Sink(
label="grid_sink"+ tag,
inputs={bel: source_flow},
)
gso = solph.Source(
label="grid_source"+ tag,
outputs={bel: solph.Flow(nominal_value=5000, variable_costs=0.1)},
)
dem = solph.Sink(
label="demand"+ tag,
inputs={bel: solph.Flow(fix=dem_data, nominal_value=10)},
)
pvi = solph.Source(
label="pv"+ tag,
outputs={bel: solph.Flow(fix=solar_data, investment=pv_invest)},
)
return (bel, gsi, gso, dem, pvi)
source_flow = solph.Flow(nominal_value=5000, variable_costs=4)
pv_invest = solph.Investment(ep_cost=0.1)
es_4 = solph.EnergySystem(timeindex=date_time_index)
es_4.add(*scenario_broken("_1", solar_data, demand_data, pv_invest, source_flow))
es_4.add(*scenario_broken("_2", solar_data, demand_data_2, pv_invest, source_flow))
So at the energy system level have flows/investments/etc… that occur at multiple buses. This doesn’t work right now (with a variety of error messages). Could it be made to work with relatively little effort? Would someone be interested in working on this with us? It seems to me like this ability would make the interface cleaner (e.g. a Python Investment object would correspond to an underlying variable) and enable a number of similar modelling approaches.
Best,
Ekaterina and Frank
The notebooks:
oemof_scenarios_minimal.ipynb (219.3 KB)
oemof_scenarios_minimal-2.ipynb (200.1 KB)
oemof_scenarios_minimal.ipynb is the system sketched above.
oemof_scenarios_minimal-2.ipynb a variant that makes the different scenarios more obvious in the solver results.