Error Message if production differs from demand

Hello everyone,
I am modeling an electrolyzer that meets demand and can offer negative balancing energy when there is capacity left.
In the model, I only look at one time period at a time and provide data from a previous optimization as input data, which determines the normal production and the balancing energy provided. I only want to check whether balancing energy is also called up.
I now have the problem that if the production does not correspond to the demand I get the error message that the problem is unsolvable or unlimited. In this case, however, the quantity produced should simply flow into the storage facility.
In this example it works if the demand is equal to 10, otherwise not.

Here is the sink where I define my demand:

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

That’s my GenericStorageComponent:

#hydrogen storage
h2_storage = solph.components.GenericStorage(
    label="hydrogen 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=storage_content_h2_storage,
    #balanced = False
    #inflow_conversion_factor=0.9,
    #outflow_conversion_factor=0.9
)

Does anyone have any idea what my mistake could be?

Here is my complete Code if needed.

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

##################### einzugebende Parameter #####################
num_tsteps = 24 #Anzahl der Zeitschritte
number_tsteps = num_tsteps

c_h2_virtual = -20 #virtueller H2-Preis €/MWh
c_h2 = 0 #angenommener H2-Preis für freien H2-Markt €/MWh
c_heat = 0
c_oxygen = 0
power_ely = 20

#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)



#sonstige Variablen
min_rate = 5/P_in_max #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.002 #loss rate per hour
el_storage_variable_costs = 0
el_storage_initial_storage_level = 0
hydrogen_storage_capacity = 150
hydrogen_storage_input_flow = 10
hydrogen_storage_output_flow = 50
hydrogen_storage_variable_costs = 0
hydrogen_storage_loss_rate = 0.002
hydrogen_storage_initial_storage_level = 0
p_demand = -10 #Preis für den H2 an feste Abnehmer verkauft wird [€/MWh]

c_el = [28.32, 10.07, -4.08, -9.91, -7.41, -12.55, -17.25, -15.07, -4.93, -6.33, -4.93, 0.45, 0.12, -0.02, 0.0, -0.03, 1.97, 9.06, 0.07, -4.97, -6.98, -24.93, -4.87, -28.93, -33.57]
lp_neg_affr = [22.2, 22.2, 22.2, 22.2, 21.52, 21.52, 21.52, 21.52, 14.84, 14.84, 14.84, 14.84, 4.04, 4.04, 4.04, 4.04, 3.12, 3.12, 3.12, 3.12, 4.8, 4.8, 4.8, 4.8, 18.24]
b_neg = [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 0, 1, 0, 1, 0, 1]

random.seed(42)
demand_h2 = [random.randint(0, 10) for _ in range(8760)] # Nachfrage

Here starts the optimization

electricity_flow_ely1_1 = 20
electricity_flow_ely1_2 = 0
electricity_flow_ely2_1 = 0
storage_content_el_storage = 0
storage_content_h2_storage = 0

%%time
n=0

c_el_angepasst = c_el.iloc[n:]#.reset_index(drop=True)
lp_neg_affr_angepasst = lp_neg_affr.iloc[n:]#.reset_index(drop=True)
demand_h2_angepasst = demand_h2[n:]



#Hier wird der Startzeitpunkt in jeder Iteration um eine Stunde nach hinten verschoben - passend zu den Listen
start_time = dt.datetime(2019, 1, 1, 0, 0) + dt.timedelta(hours=n)

#definition of time index
datetime_index = solph.create_time_index(year=2019, number=1, start=start_time)

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

#Definition Bus-Components
b_el = solph.Bus("electricity bus")
b_h2 = solph.Bus("hydrogen bus")
b_heat = solph.Bus("heat bus")
b_o2 = solph.Bus("oxygen bus")
b_el_neg_affr = solph.Bus("neg affr electricity 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_angepasst
        )
    }
)



source_el_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_angepasst,
            nonconvex=solph.NonConvex(
            )            
        )
    }
)

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

#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
        )
    }
)


#### Electrolyzer hydrogen market ####
#firt part electrolyzer to cover hydrogen demand/production
electrolyzer1_1 = solph.components.OffsetConverter(
    label='electrolyzer market 1',
    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
    }
)



#fourth part electrolyzer to cover delivering of neg. balancing energy
electrolyzer2_2 = 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_h2: solph.Flow()
    },
    conversion_factors={
        b_h2: slope_h2,
    },
    normed_offsets={
        b_h2: offset_h2,
    }
)



#### Storages ####
#battery storage
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=storage_content_el_storage,
    #balanced = False
    #inflow_conversion_factor=0.9,
    #outflow_conversion_factor=0.9
)

#hydrogen storage
h2_storage = solph.components.GenericStorage(
    label="hydrogen 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=storage_content_h2_storage,
    #balanced = False
    #inflow_conversion_factor=0.9,
    #outflow_conversion_factor=0.9
)


es3.add(b_el, b_h2, b_heat, b_o2, b_el_neg_affr,
        source_el, sink_h2_demand, sink_heat, sink_o2, source_el_neg_affr,
        electrolyzer1_1, el_storage, h2_storage, electrolyzer2_2
)

om3 = solph.Model(es3)

myblock3 = po.Block()
om3.add_component("MyBlock3", myblock3)

def flow_ely1_1(m,t):
    expr = om3.flow[b_el, electrolyzer1_1, t] == electricity_flow_ely1_1
    return expr

myblock3.flow_ely1_1 = po.Constraint(om3.TIMESTEPS, rule=flow_ely1_1)



#at call-off of neg. balacing energy (b_neg=1) flow of ely2_1 and ely2_2 have to be equal
def production_constraint_neg_affr_1(m,t):
    expr = b_neg[n] * om3.flow[b_el_neg_affr, electrolyzer2_2, t] == b_neg[n] * electricity_flow_ely2_1
    return expr

myblock3.production_constraint_neg_affr_1 = po.Constraint(om3.TIMESTEPS, rule=production_constraint_neg_affr_1)


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

myblock3.production_constraint_neg_affr_2 = po.Constraint(om3.TIMESTEPS, rule=production_constraint_neg_affr_2)



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

results3 = solph.views.convert_keys_to_strings(om3.results(), keep_none_type=True)

Thanks for any idea :slight_smile: