Planning delay of a flow

[dhorschi]

Hello everyone,

I am modeling an electrolyzer that participates in the balancing market.
I am creating several transformers to model the electrolyzer on the different markets.
Now I would like to achieve that the hydrogen which is additionally produced by negative balancing energy is always available one period after the actual production.
My idea for this was to insert an additional storage in which the produced hydrogen has to remain for at least one period of time.

I saw the parameter lifetime_outflow in the docs, but I don’t think it quite fits my problem.

Does anyone have an idea how I can implement my idea? Other ideas for my problem are welcome as well.

Thanks for your ideas
Kind regards

Luca

[dhorschi]

Hi everyone,

I think I have found a way and wanted to share it with you.

t_delay = 1
def delay_constraint(m,t):
    if t==0:
        expr = om.flow[help_storage, b_h2, t] == 0
        return expr
    elif t >= t_delay:
        expr = om.flow[help_storage, b_h2, t] == om.flow[b_h2_help, help_storage, t-t_delay]
        return expr
    else:
        return po.Constraint.Skip

myblock.delay_constraint = po.Constraint(om.TIMESTEPS, rule=delay_constraint)

Unfortunately, this brings the next problem
In my model, I enter a list with zeros and ones to indicate whether control energy is called up or not. Thats why the optimizer already knows whether hydrogen would be produced by calling up balancing energy in every period. So I want to delay the use of this hydrogen by one period, since in reality the demand would have to be covered independently of the balancing energy. To do this, I use a “help storage” and the corresponding constraint.

With the constraint, however, no more balancing energy is provided and I don’t understand why. Even if I set the capacity price unrealistically high.

I’m sharing my code so that you can hopefully understand the problem better. Unfortunately, this makes the message quite long.

from oemof import solph
import pandas as pd
import numpy as np
import pyomo.environ as po
import random

##################### einzugebende Parameter #####################
num_tsteps = 8 #Anzahl der Zeitschritte
c_h2_virtual = -20 #virtueller H2-Preis €/MWh
c_h2 = 1000 #angenommener H2-Preis für freien H2-Markt €/MWh
c_heat = 0
c_oxygen = 0

power_ely = 20  # Beispiel-Leistung des Elektrolyseurs in MW
investment_cost_ely = 1500  # Beispiel-Investitionskosten in Euro/kW
number_years = 20 #Abschreibungsdauer in Jahren
interest_rate = 0.06 #Zinssatz für die Abschreibung

#annualized capex
a = (((1+interest_rate)**number_years)*interest_rate)/(((1+interest_rate)**number_years)-1)
annualized_cost = investment_cost_ely * power_ely * 1000 * a


#Daten Teillastverahlten
eta_h2_min = 0.60 #efficiency at P_in_min
eta_h2_max = 0.50 #efficiency at P_in_max
P_in_min = 2 
P_in_max = 20

#slope und offset for part-load behavior hydrogen, heat, oxygen
slope_h2, offset_h2 = solph.components._offset_converter.slope_offset_from_nonconvex_input(P_in_max, P_in_min, eta_at_max=0.5, eta_at_min=0.6)

slope_heat, offset_heat = solph.components._offset_converter.slope_offset_from_nonconvex_input(P_in_max, P_in_min, eta_at_max=0.4, eta_at_min=0.3)

slope_o2, offset_o2 = solph.components._offset_converter.slope_offset_from_nonconvex_input(P_in_max, P_in_min, eta_at_max=0.2, eta_at_min=0.1)

#Daten Teillastverahlten
eta_h2_min_2 = 0 #efficiency at P_in_min
eta_h2_max_2 = 0.60 #efficiency at P_in_max
P_in_min_2 = 0 
P_in_max_2 = 2

#slope und offset for part-load behavior hydrogen, heat, oxygen
slope_h2_2, offset_h2_2 = solph.components._offset_converter.slope_offset_from_nonconvex_input(P_in_max_2, P_in_min_2, eta_at_max=eta_h2_max_2, eta_at_min=eta_h2_min_2)

slope_heat_2, offset_heat_2 = solph.components._offset_converter.slope_offset_from_nonconvex_input(P_in_max_2, P_in_min_2, eta_at_max=0.4, eta_at_min=0)

slope_o2_2, offset_o2_2 = solph.components._offset_converter.slope_offset_from_nonconvex_input(P_in_max_2, P_in_min_2, eta_at_max=0.1, eta_at_min=0)


#sonstige Variablen
min_rate = 5/power_ely #Rate damit minimale Regelenergie 5MW ist
el_storage_capacity = 5
el_storage_input_flow = 5
el_storage_output_flow = 5
el_storage_loss_rate = 0.02 #loss rate per hour
el_storage_variable_costs = 0
el_storage_initial_storage_level = 0
hydrogen_storage_capacity = 50
hydrogen_storage_input_flow = 5
hydrogen_storage_output_flow = 5
hydrogen_storage_variable_costs = 0
hydrogen_storage_loss_rate = 0.02
hydrogen_storage_initial_storage_rate = 0
p_demand = -10 #Preis für den H2 an feste Abnehmer verkauft wird [€/MWh]

c_el=20 # electricity price in €/MWh
b_neg = [0,1,1,1,1,1,1,1] #if value is 1 - balancing energy is called up
lp_neg_affr = 200 # capacity price in €/MW 
ap_neg_affr= 30 #working price in €/MWh

demand_h2 = [6,1,4,9,8,3,4,5]


#definition of time index
datetime_index = solph.create_time_index(year=2023, number=num_tsteps)

es2 = solph.EnergySystem(timeindex=datetime_index, infer_last_interval=False)

#Definition Bus-Components
b_el = solph.Bus("electricity bus")
b_el_neg_affr = solph.Bus("neg affr electricity")
b_el_free = solph.Bus("free electricity bus")
b_h2 = solph.Bus("hydrogen bus")
b_h2_neg_affr_virt = solph.Bus("neg affr hydrogen bus virt")
b_heat = solph.Bus("heat bus")
b_o2 = solph.Bus("oxygen bus")
b_h2_help = solph.Bus("hydrogen help bus")

##### Definition der Komponenten #####

#electricity source for basic hydrogen demand
source_el = solph.components.Source(
    "electricity import", 
    outputs={
        b_el: solph.Flow(
            variable_costs=c_el
        )
    }
)

#second electricity source for call-off electrolyzer for neg. balancing energy
source_neg_affr = solph.components.Source(
    "electricity neg affr",
    outputs={
        b_el_neg_affr: solph.Flow(
            nominal_value=power_ely, 
            min=min_rate, 
            variable_costs=c_el,
            nonconvex=solph.NonConvex()
            #nonconvex=solph.NonConvex(
                #minimum_downtime=4, initial_status=0 
            #)
        )
    }
)

#third electricity source for keepin neg. balancing energy available
source_el_free = solph.components.Source(
    "free electricity",
    outputs={
        b_el_free: solph.Flow(
            nominal_value=power_ely,
            min=min_rate,
            nonconvex=solph.NonConvex(
                #minimum_downtime=4, #initial_status=0
            )            
        )
    }
)

#Sink for fix haydrogen demand via contract
sink_h2_demand = solph.components.Sink(
    "hydrogen demand",
    inputs={
        b_h2: solph.Flow(
            fix=demand_h2,
            nominal_value=1,
            variable_costs=c_h2_virtual
        )
    }
)

#Sink for selling hydrogen at a potential market
sink_h2_market = solph.components.Sink(
    "hydrogen market",
    inputs={
        b_h2: solph.Flow(
            variable_costs=c_h2 
        )
    }
)

sink_h2_neg_affr_virt = solph.components.Sink(
    "neg affr h2 sink virt",
    inputs={
        b_h2_neg_affr_virt: solph.Flow(
            variable_costs= -lp_neg_affr #price for keeping neg. balancing energy available 
        )        
    }
)

#sink for byproduct heat
sink_heat = solph.components.Sink(
    "heat export",
    inputs={
        b_heat: solph.Flow(
            variable_costs=c_heat
        )
    }
)


#sink for byproduct oxygen
sink_o2 = solph.components.Sink(
    "oxygen export",
    inputs={
        b_o2: solph.Flow(
            variable_costs=c_oxygen
        )
    }
)


#firt part electrolyzer to cover hydrogen demand/production
electrolyzer1 = solph.components.OffsetConverter(
    label='electrolyzer market',
    inputs={
        b_el: solph.Flow(
            nominal_value=P_in_max,
            nonconvex=solph.NonConvex(),
            min=P_in_min / P_in_max,
        )
    },
    outputs={
        b_heat: solph.Flow(),
        b_h2: solph.Flow(),
        b_o2: solph.Flow(),
    },
    conversion_factors={
        b_heat: slope_heat,
        b_h2: slope_h2,
        b_o2: slope_o2
    },
    normed_offsets={
        b_heat: offset_heat,
        b_h2: offset_h2,
        b_o2: offset_o2
    }
)

#firt part electrolyzer to cover hydrogen demand/production
electrolyzer1_2 = solph.components.OffsetConverter(
    label='electrolyzer market 2',
    inputs={
        b_el: solph.Flow(
            nominal_value=P_in_max_2,
            nonconvex=solph.NonConvex(),
            min=P_in_min_2 / P_in_max_2,
        )
    },
    outputs={
        b_heat: solph.Flow(),
        b_h2: solph.Flow(),
        b_o2: solph.Flow(),
    },
    conversion_factors={
        b_heat: slope_heat_2,
        b_h2: slope_h2_2,
        b_o2: slope_o2_2
    },
    normed_offsets={
        b_heat: offset_heat_2,
        b_h2: offset_h2_2,
        b_o2: offset_o2_2
    }
)



#second part electrolyzer to cover holding of neg. balancing energy
electrolyzer2 = solph.components.OffsetConverter(
    label='electrolyzer neg affr holding',
    inputs={
        b_el_free: solph.Flow(
            nominal_value=P_in_max,
            nonconvex=solph.NonConvex(),
            min=P_in_min / P_in_max,
        )
    },
    outputs={
        b_h2_neg_affr_virt: solph.Flow()
    },
    conversion_factors={
        b_h2_neg_affr_virt: 1 #slope_h2,
    },
    normed_offsets={
        b_h2_neg_affr_virt: 0 #offset_h2,
    }
)

#third part electrolyzer to cover call-off of neg. balancing energy
electrolyzer2_1 = solph.components.OffsetConverter(
    label='electrolyzer neg affr delivering',
    inputs={
        b_el_neg_affr: solph.Flow(
            nominal_value=P_in_max,
            nonconvex=solph.NonConvex(),
            min=P_in_min / P_in_max,
        )
    },
    outputs={
        b_heat: solph.Flow(),
        #b_h2_neg_affr: solph.Flow(
        b_h2_help:solph.Flow(    
            variable_costs=-ap_neg_affr
        ),
        b_o2: solph.Flow(),
    },
    conversion_factors={
        b_heat: slope_heat,
        b_h2_help: slope_h2,
        #b_h2_neg_affr: slope_h2,
        b_o2: slope_o2        
    },
    normed_offsets={
        b_heat: offset_heat,
        b_h2_help: offset_h2,
        #b_h2_neg_affr: offset_h2,
        b_o2: offset_o2,
    }
)

el_storage = solph.components.GenericStorage(
    label="electricity storage",
    nominal_storage_capacity = el_storage_capacity,
    inputs={
        b_el: solph.Flow(
            nominal_value=el_storage_input_flow,
            variable_costs=el_storage_variable_costs,
            nonconvex=solph.NonConvex()
        )
    },
    outputs={
        b_el:solph.Flow(
            nominal_value=el_storage_output_flow
        )
    },
    loss_rate=el_storage_loss_rate,
    initial_storage_level=el_storage_initial_storage_level,
    #inflow_conversion_factor=0.9,
    #outflow_conversion_factor=0.9
)


storage = solph.components.GenericStorage(
    label="storage",
    nominal_storage_capacity = hydrogen_storage_capacity,
    inputs={
        b_h2: solph.Flow(
            nominal_value=hydrogen_storage_input_flow,
            variable_costs=hydrogen_storage_variable_costs,
            nonconvex=solph.NonConvex()
        )
    },
    outputs={
        b_h2:solph.Flow(
            nominal_value=hydrogen_storage_output_flow
        )
    },
    loss_rate=hydrogen_storage_loss_rate,
    initial_storage_level=hydrogen_storage_initial_storage_rate,
    #inflow_conversion_factor=0.9,
    #outflow_conversion_factor=0.9
)

help_storage = solph.components.GenericStorage(
    label="help storage",
    nominal_storage_capacity = hydrogen_storage_capacity,
    inputs={
        b_h2_help: solph.Flow(
            nominal_value=hydrogen_storage_input_flow,
            variable_costs=hydrogen_storage_variable_costs,
            nonconvex=solph.NonConvex()
        )
    },
    outputs={
        b_h2:solph.Flow(
            nominal_value=hydrogen_storage_output_flow
        )
    },
    loss_rate=hydrogen_storage_loss_rate,
    initial_storage_level=hydrogen_storage_initial_storage_rate,
    #inflow_conversion_factor=0.9,
    #outflow_conversion_factor=0.9
)



es2.add(
    b_el, b_heat, b_h2, b_o2, b_el_neg_affr,  b_el_free, b_h2_neg_affr_virt, b_h2_help, #b_el_pos_affr, b_h2_pos_affr,
    source_el, sink_h2_demand, sink_h2_market, sink_heat, sink_o2, source_neg_affr, sink_h2_neg_affr_virt, source_el_free, #source_pos_affr, sink_pos_affr, 
    electrolyzer1, electrolyzer2, electrolyzer2_1,  el_storage, storage, electrolyzer1_2, help_storage #electrolyzer3,
)

om = solph.Model(es2)



                            ########################### adding Custom Constraint ########################### 

myblock = po.Block()
om.add_component("MyBlock", myblock)

solph.constraints.limit_active_flow_count(
        om, "flow count", [(b_el, electrolyzer1), (b_el, electrolyzer1_2)], lower_limit=0, upper_limit=1
    )


#limit the maximung power of the part electrolyzers
def limit_active_flow_count_rule(m, t):
    expr = (om.flow[b_el, electrolyzer1, t] + om.flow[b_el_free, electrolyzer2, t] + om.flow[b_el, electrolyzer1_2, t]  <= P_in_max) 
    return expr

myblock.limit_active_flow_count = po.Constraint(om.TIMESTEPS, rule=limit_active_flow_count_rule)

#at call-off of neg. balacing energy (b_neg=1) flow of ely2 and ely2_1 have to be equal
def production_constraint_neg_affr_1(m,t):
    expr = b_neg[t] * om.flow[b_el_neg_affr, electrolyzer2_1, t] == b_neg[t] * om.flow[b_el_free, electrolyzer2, t]
    return expr

myblock.production_constraint_neg_affr_1 = po.Constraint(om.TIMESTEPS, rule=production_constraint_neg_affr_1)


#if there is no call-off of neg. balancing energy (b_neg=0)the production of ely2_1 has to be 0 - DOES NOT WORK
def production_constraint_neg_affr_2(m,t):
    expr = (1-b_neg[t]) * om.flow[b_el_neg_affr, electrolyzer2_1, t] == 0
    return expr

myblock.production_constraint_neg_affr_2 = po.Constraint(om.TIMESTEPS, rule=production_constraint_neg_affr_2)

#DateTimeIndex with start times of the control energy time slices
decision_times = pd.date_range(datetime_index[0], datetime_index[-1], freq='4h')

def time_constraint_neg_affr(m,t):
    if datetime_index[t] not in decision_times:
        expr = om.flow[source_el_free, b_el_free, t] == om.flow[source_el_free, b_el_free,t-1]
        return expr
    else:
        return po.Constraint.Skip

myblock.time_constraint_neg_affr = po.Constraint(om.TIMESTEPS, rule=time_constraint_neg_affr)


t_delay = 1
def delay_constraint(m,t):
    if t==0:
        expr = om.flow[help_storage, b_h2, t] == 0
        return expr
    elif t >= t_delay:
        expr = om.flow[help_storage, b_h2, t] == om.flow[b_h2_help, help_storage, t-t_delay]
        return expr
    else:
        return po.Constraint.Skip

myblock.delay_constraint = po.Constraint(om.TIMESTEPS, rule=delay_constraint)


#lösen des Optimierungsproblems
om.solve("gurobi")

results = solph.views.convert_keys_to_strings(om.results(), keep_none_type=True)

It would be great if someone could check my code. Maybe @pschoen you have an idea since you’re helping so many in here.

Thank you very much in advance

[pschoen]

Hi @dhorschi,

since you mention me: Many questions asked here, are very similar to things I at least thought about. So, I can answer right after reading. This one is interesting, but new to me. Also, it’s to complex for me to completely get all the consequences at one glance.

What I see, however, is that your help_storage has a non-zero loss rate. At the same time, your constraint defines that inflow and outflow must be the same. As a consequence, the flow is forced to be zero: If there was a non-zero inflow, there would be losses, so the outflow would have to be smaller.

I would suggest that you divide and conquer your problem. The constraint itself should do what it is supposed to do, but maybe check it. (I am really bad at this, I always need to debug stuff.) To simplify your constrained problem, you can replace your storage by a slack bus (Bus(balanced=False)) between a source and a sink. Now, if you force-push (Flow(nominal_value=, fix=[...])) something into the central slack bus, it should come out to the sink one step later due to the constraint. If that works, you can add that to your more complex model.

Cheers,
Patrik