diff --git a/src/topupopt/problems/esipp/blocks/delta.py b/src/topupopt/problems/esipp/blocks/delta.py
new file mode 100644
index 0000000000000000000000000000000000000000..59b418feb43902d31282c8666f7417ce3071b7f1
--- /dev/null
+++ b/src/topupopt/problems/esipp/blocks/delta.py
@@ -0,0 +1,369 @@
+# imports
+
+import pyomo.environ as pyo
+    
+# *****************************************************************************
+# *****************************************************************************
+
+def price_delta_block(
+    model: pyo.AbstractModel,
+    enable_default_values: bool = True,
+    enable_validation: bool = True,
+    enable_initialisation: bool = True
+    ):
+    
+    # auxiliary set for pyomo
+    model.set_GLQPK = model.set_GL_exp_imp*model.set_QPK
+    
+    # create a block for a transshipment node during a given time interval
+    def rule_block_prices(b, g, l, q, p, k):
+        
+        # *********************************************************************
+        # *********************************************************************
+        
+        # sets
+        
+        # set of price segments
+        b.set_S = pyo.Set()
+        
+        # TODO: introduce a set of price segments for non-convex tariffs
+        # def init_set_S_nonconvex(b, s):
+        #     return (s for s in b.set_S if b.param_price_function_is_convex)
+        # b.set_S_nonconvex = pyo.Set(within=b.set_S, initialize=init_set_S_nonconvex)
+    
+        # *********************************************************************
+        # *********************************************************************
+        
+        # parameters
+
+        # resource prices
+        b.param_p_s = pyo.Param(b.set_S, within=pyo.NonNegativeReals)
+        
+        # price function convexity
+        b.param_price_function_is_convex = pyo.Param(within=pyo.Boolean)
+
+        # maximum resource volumes for each prices
+        b.param_v_max_s = pyo.Param(b.set_S, within=pyo.NonNegativeReals)
+    
+        # *********************************************************************
+        # *********************************************************************
+    
+        # variables
+    
+        # *********************************************************************
+        # *********************************************************************
+        
+        # TODO: consider replacing None in the bounds with inf in the parameter
+    
+        # import and export flows
+        def bounds_var_trans_flows(b):
+            try:
+                return (0, sum(b.param_v_max_s[s] for s in b.set_S))
+            except Exception:
+                return (0, None)
+        b.var_trans_flows = pyo.Var(
+            # within=pyo.NonNegativeReals,
+            bounds=bounds_var_trans_flows
+        )
+        
+        # import and export flows
+        def rule_var_segment_usage_s(b, s):
+            if b.param_price_function_is_convex:
+                return pyo.UnitInterval
+            else:
+                return pyo.Binary
+        b.var_segment_usage_s = pyo.Var(
+            b.set_S, 
+            within=rule_var_segment_usage_s, 
+        )
+        
+        # *********************************************************************
+        # *********************************************************************
+    
+        # import flow costs and export flow revenues
+        def rule_constr_trans_monetary_flows(b):
+            if (g,l) in b.parent_block().set_GL_imp:
+                return (
+                    sum(
+                        (b.var_segment_usage_s[s]*b.param_p_s[s]*b.param_v_max_s[s] 
+                         if s == 0 else
+                         b.var_segment_usage_s[s]*(b.param_p_s[s]*b.param_v_max_s[s]-b.param_p_s[s-1]*b.param_v_max_s[s-1]))
+                        for s in b.set_S
+                    )
+                    == b.parent_block().var_ifc_glqpk[(g,l,q,p,k)]
+                )
+            else:
+                return (
+                    sum(
+                        (b.var_segment_usage_s[s]*b.param_p_s[s]*b.param_v_max_s[s] 
+                         if s == 0 else
+                         b.var_segment_usage_s[s]*(b.param_p_s[s]*b.param_v_max_s[s]-b.param_p_s[s-1]*b.param_v_max_s[s-1]))
+                        for s in b.set_S
+                    )
+                    == b.parent_block().var_efr_glqpk[(g,l,q,p,k)]
+                )
+        b.constr_trans_monetary_flows = pyo.Constraint(
+            rule=rule_constr_trans_monetary_flows
+        )        
+        
+        # imported and exported flows: volumes
+        def rule_constr_trans_flows_seg(b):
+            return sum(
+                b.var_segment_usage_s[s]*b.param_v_max_s[s]
+                for s in b.set_S
+            ) == b.var_trans_flows
+        b.constr_trans_flows_seg = pyo.Constraint(rule=rule_constr_trans_flows_seg)
+        
+        # imported and exported flows: energy system
+        def rule_constr_trans_flows_sys(b):
+            if (g,l) in b.parent_block().set_GL_imp:
+                return sum(
+                    b.parent_block().var_v_glljqk[(g,l,l_star,j,q,k)]
+                    for l_star in b.parent_block().set_L[g]
+                    if l_star not in b.parent_block().set_L_imp[g]
+                    for j in b.parent_block().set_J[(g,l,l_star)]  # only directed arcs
+                ) == b.var_trans_flows
+            else:
+                return sum(
+                    b.parent_block().var_v_glljqk[(g,l_star,l,j,q,k)]
+                    * b.parent_block().param_eta_glljqk[(g,l_star,l,j,q,k)]
+                    for l_star in b.parent_block().set_L[g]
+                    if l_star not in b.parent_block().set_L_exp[g]
+                    for j in b.parent_block().set_J[(g,l_star,l)]  # only directed arcs
+                ) == b.var_trans_flows
+        b.constr_trans_flows_sys = pyo.Constraint(rule=rule_constr_trans_flows_sys)
+    
+        # *********************************************************************
+        # *********************************************************************
+        
+        # non-convex price functions
+        if not b.param_price_function_is_convex:
+            # this means the segment usage variables are binary
+            
+            raise NotImplementedError
+            
+            # # if delta var is one, previous ones must be one too
+            # # if delta var is zero, the next ones must also be zero
+            # def rule_constr_delta_summing_logic(b, s):
+            #     if s == len(b.set_S)-1:
+            #         # last segment, skip
+            #         # convex, skip
+            #         return pyo.Constraint.Skip
+            #     return (
+            #         b.var_segment_usage_s[s] >= 
+            #         b.var_segment_usage_s[s+1]
+            #         )
+            # b.constr_delta_summing_logic = pyo.Constraint(
+            #     b.set_S, rule=rule_constr_delta_summing_logic
+            #     )
+        
+        # *********************************************************************
+        # ********************************************************************* 
+        
+    model.block_prices = pyo.Block(model.set_GLQPK, rule=rule_block_prices)
+
+    # *************************************************************************
+    # *************************************************************************
+
+# *****************************************************************************
+# *****************************************************************************
+# *****************************************************************************
+# *****************************************************************************
+# *****************************************************************************
+# *****************************************************************************
+# *****************************************************************************
+# *****************************************************************************
+# *****************************************************************************
+# *****************************************************************************
+# *****************************************************************************
+# *****************************************************************************
+
+def price_delta_no_block(
+    model: pyo.AbstractModel,
+    # convex_price_function: bool = True,
+    enable_default_values: bool = True,
+    enable_validation: bool = True,
+    enable_initialisation: bool = True
+    ):
+    
+    raise NotImplementedError
+    
+    # auxiliary set for pyomo
+    model.set_GLQPK = model.set_GL_exp_imp*model.set_QPK
+
+    # set of price segments
+    model.set_S = pyo.Set(model.set_GLQPK)
+
+    # set of GLQKS tuples
+    def init_set_GLQPKS(m):
+        return (
+            (g, l, q, p, k, s)
+            # for (g,l) in m.set_GL_exp_imp
+            # for (q,k) in m.set_QK
+            for (g, l, q, p, k) in m.set_S
+            for s in m.set_S[(g, l, q, p, k)]
+        )
+
+    model.set_GLQPKS = pyo.Set(
+        dimen=6, initialize=(init_set_GLQPKS if enable_initialisation else None)
+    )
+
+    # *************************************************************************
+    # *************************************************************************
+
+    # parameters
+
+    # resource prices
+
+    model.param_p_glqpks = pyo.Param(model.set_GLQPKS, within=pyo.NonNegativeReals)
+    
+    # price function convexity
+
+    model.param_price_function_is_convex = pyo.Param(
+        model.set_GLQPK, 
+        within=pyo.Boolean
+        )
+
+    # maximum resource volumes for each prices
+
+    model.param_v_max_glqpks = pyo.Param(
+        model.set_GLQPKS, 
+        within=pyo.NonNegativeReals,
+        # default=math.inf
+        )
+
+    # *************************************************************************
+    # *************************************************************************
+
+    # variables
+
+    # *************************************************************************
+    # *************************************************************************
+
+    # import and export flows
+    def bounds_var_trans_flows_glqpks(m, g, l, q, p, k, s):
+        # return (0, m.param_v_max_glqpks[(g, l, q, p, k, s)])
+        if (g, l, q, p, k, s) in m.param_v_max_glqpks:
+            # predefined finite capacity
+            return (0, m.param_v_max_glqpks[(g, l, q, p, k, s)])
+        else:
+            # infinite capacity
+            return (0, None)
+    model.var_trans_flows_glqpks = pyo.Var(
+        model.set_GLQPKS, within=pyo.NonNegativeReals, bounds=bounds_var_trans_flows_glqpks
+    )
+
+    # *************************************************************************
+    # *************************************************************************
+
+    # import flow costs and export flow revenues
+    def rule_constr_trans_monetary_flows(m, g, l, q, p, k):
+        if (g,l) in m.set_GL_imp:
+            return (
+                sum(
+                    m.var_trans_flows_glqpks[(g, l, q, p, k, s)]
+                    * m.param_p_glqpks[(g, l, q, p, k, s)]
+                    for s in m.set_S[(g, l, q, p, k)]
+                )
+                == m.var_ifc_glqpk[(g, l, q, p, k)]
+            )
+        else:
+            return (
+                sum(
+                    m.var_trans_flows_glqpks[(g, l, q, p, k, s)]
+                    * m.param_p_glqpks[(g, l, q, p, k, s)]
+                    for s in m.set_S[(g, l, q, p, k)]
+                )
+                == m.var_efr_glqpk[(g, l, q, p, k)]
+            )
+    model.constr_trans_monetary_flows = pyo.Constraint(
+        model.set_GLQPK, rule=rule_constr_trans_monetary_flows
+    )
+
+    # imported and exported flows
+    def rule_constr_trans_flows(m, g, l, q, p, k):
+        if (g,l) in m.set_GL_imp:
+            return sum(
+                m.var_v_glljqk[(g, l, l_star, j, q, k)]
+                for l_star in m.set_L[g]
+                if l_star not in m.set_L_imp[g]
+                for j in m.set_J[(g, l, l_star)]  # only directed arcs
+            ) == sum(m.var_trans_flows_glqpks[(g, l, q, p, k, s)] for s in m.set_S[(g, l, q, p, k)])
+        else:
+            return sum(
+                m.var_v_glljqk[(g, l_star, l, j, q, k)]
+                * m.param_eta_glljqk[(g, l_star, l, j, q, k)]
+                for l_star in m.set_L[g]
+                if l_star not in m.set_L_exp[g]
+                for j in m.set_J[(g, l_star, l)]  # only directed arcs
+            ) == sum(m.var_trans_flows_glqpks[(g, l, q, p, k, s)] for s in m.set_S[(g, l, q, p, k)])
+
+    model.constr_trans_flows = pyo.Constraint(
+        model.set_GLQPK, rule=rule_constr_trans_flows
+    )
+
+    # *************************************************************************
+    # *************************************************************************
+    
+    # non-convex price functions
+    # TODO: remove these variables from the model if they are not needed
+    # delta variables
+    model.var_active_segment_glqpks = pyo.Var(
+        model.set_GLQPKS, within=pyo.Binary
+    )
+    
+    # segments must be empty if the respective delta variable is zero
+    def rule_constr_empty_segment_if_delta_zero(m, g, l, q, p, k, s):
+        if len(m.set_S[(g,l,q,p,k)]) == 1 or m.param_price_function_is_convex[(g,l,q,p,k)]:
+            # single segment, skip
+            # convex, skip
+            return pyo.Constraint.Skip
+        return (
+            m.var_trans_flows_glqpks[(g,l,q,p,k,s)] <= 
+            m.param_v_max_glqpks[(g,l,q,p,k,s)]*
+            m.var_active_segment_glqpks[(g,l,q,p,k,s)]
+            )
+    model.constr_empty_segment_if_delta_zero = pyo.Constraint(
+        model.set_GLQPKS, rule=rule_constr_empty_segment_if_delta_zero
+        )
+    
+    # if delta var is one, previous ones must be one too
+    # if delta var is zero, the next ones must also be zero
+    def rule_constr_delta_summing_logic(m, g, l, q, p, k, s):
+        if s == len(m.set_S[(g,l,q,p,k)])-1 or m.param_price_function_is_convex[(g,l,q,p,k)]:
+            # last segment, skip
+            # convex, skip
+            return pyo.Constraint.Skip
+        return (
+            m.var_active_segment_glqpks[(g,l,q,p,k,s)] >= 
+            m.var_active_segment_glqpks[(g,l,q,p,k,s+1)]
+            )
+    model.constr_delta_summing_logic = pyo.Constraint(
+        model.set_GLQPKS, rule=rule_constr_delta_summing_logic
+        )
+    
+    # if a segment is not completely used, the next ones must remain empty
+    def rule_constr_fill_up_segment_before_next(m, g, l, q, p, k, s):
+        if s == len(m.set_S[(g,l,q,p,k)])-1 or m.param_price_function_is_convex[(g,l,q,p,k)]:
+            # last segment, skip
+            # convex, skip
+            return pyo.Constraint.Skip
+        return (
+            m.var_trans_flows_glqpks[(g,l,q,p,k,s)] >= 
+            m.var_active_segment_glqpks[(g,l,q,p,k,s+1)]*
+            m.param_v_max_glqpks[(g,l,q,p,k,s)]
+            )
+        # return (
+        #     m.var_if_glqpks[(g,l,q,p,k,s)]/m.param_v_max_glqpks[(g,l,q,p,k,s)] >= 
+        #     m.var_active_segment_glqpks[(g,l,q,p,k,s+1)]
+        #     )
+        # return (
+        #     m.param_v_max_glqpks[(g,l,q,p,k,s)]-m.var_if_glqpks[(g,l,q,p,k,s)] <= 
+        #     m.param_v_max_glqpks[(g,l,q,p,k,s)]*(1- m.var_active_segment_glqpks[(g,l,q,p,k,s+1)])
+        #     )
+    model.constr_fill_up_segment_before_next = pyo.Constraint(
+        model.set_GLQPKS, rule=rule_constr_fill_up_segment_before_next
+        )
+
+# *****************************************************************************
+# *****************************************************************************
\ No newline at end of file
diff --git a/src/topupopt/problems/esipp/blocks/lambda.py b/src/topupopt/problems/esipp/blocks/lambda.py
new file mode 100644
index 0000000000000000000000000000000000000000..c09943f583629ecae62edc3382a163b0d681b7a9
--- /dev/null
+++ b/src/topupopt/problems/esipp/blocks/lambda.py
@@ -0,0 +1,388 @@
+# imports
+
+import pyomo.environ as pyo
+    
+# *****************************************************************************
+# *****************************************************************************
+
+def price_lambda_block(
+    model: pyo.AbstractModel,
+    enable_default_values: bool = True,
+    enable_validation: bool = True,
+    enable_initialisation: bool = True
+    ):
+    
+    # auxiliary set for pyomo
+    model.set_GLQPK = model.set_GL_exp_imp*model.set_QPK
+    
+    # create a block for a transshipment node during a given time interval
+    def rule_block_prices(b, g, l, q, p, k):
+        
+        # *********************************************************************
+        # *********************************************************************
+        
+        # sets
+        
+        # set of price segments
+        b.set_S = pyo.Set()
+        
+        # TODO: introduce a set of price segments for non-convex tariffs
+        # def init_set_S_nonconvex(b, s):
+        #     return (s for s in b.set_S if b.param_price_function_is_convex)
+        # b.set_S_nonconvex = pyo.Set(within=b.set_S, initialize=init_set_S_nonconvex)
+    
+        # *********************************************************************
+        # *********************************************************************
+        
+        # parameters
+
+        # resource prices
+        b.param_p_s = pyo.Param(b.set_S, within=pyo.NonNegativeReals)
+        
+        # price function convexity
+        b.param_price_function_is_convex = pyo.Param(within=pyo.Boolean)
+
+        # maximum resource volumes for each prices
+        b.param_v_max_s = pyo.Param(b.set_S, within=pyo.NonNegativeReals)
+    
+        # *********************************************************************
+        # *********************************************************************
+    
+        # variables
+    
+        # *********************************************************************
+        # *********************************************************************
+        
+        # TODO: consider replacing None in the bounds with inf in the parameter
+    
+        # import and export flows
+        def bounds_var_trans_flows_s(b, s):
+            if s in b.param_v_max_s:
+                # predefined finite capacity
+                return (0, b.param_v_max_s[s])
+            else:
+                # infinite capacity
+                return (0, None)
+        b.var_trans_flows_s = pyo.Var(
+            b.set_S, 
+            within=pyo.NonNegativeReals, 
+            bounds=bounds_var_trans_flows_s
+        )
+    
+        # import and export flows
+        def bounds_var_trans_flows_s(b, s):
+            if s in b.param_v_max_s:
+                # predefined finite capacity
+                return (0, b.param_v_max_s[s])
+            else:
+                # infinite capacity
+                return (0, None)
+        b.var_trans_flows_s = pyo.Var(
+            b.set_S, 
+            within=pyo.NonNegativeReals, 
+            bounds=bounds_var_trans_flows_s
+        )
+        
+        # *********************************************************************
+        # *********************************************************************
+    
+        # import flow costs and export flow revenues
+        def rule_constr_trans_monetary_flows(b):
+            if (g,l) in b.parent_block().set_GL_imp:
+                return (
+                    sum(b.var_trans_flows_s[s]*b.param_p_s[s]
+                        for s in b.set_S)
+                    == b.parent_block().var_ifc_glqpk[(g,l,q,p,k)]
+                )
+            else:
+                return (
+                    sum(b.var_trans_flows_s[s]*b.param_p_s[s]
+                        for s in b.set_S)
+                    == b.parent_block().var_efr_glqpk[(g,l,q,p,k)]
+                )
+        b.constr_trans_monetary_flows = pyo.Constraint(
+            rule=rule_constr_trans_monetary_flows
+        )
+        
+        # imported and exported flows
+        def rule_constr_trans_flows(b):
+            if (g,l) in b.parent_block().set_GL_imp:
+                return sum(
+                    b.parent_block().var_v_glljqk[(g,l,l_star,j,q,k)]
+                    for l_star in b.parent_block().set_L[g]
+                    if l_star not in b.parent_block().set_L_imp[g]
+                    for j in b.parent_block().set_J[(g,l,l_star)]  # only directed arcs
+                ) == sum(b.var_trans_flows_s[s] for s in b.set_S)
+            else:
+                return sum(
+                    b.parent_block().var_v_glljqk[(g,l_star,l,j,q,k)]
+                    * b.parent_block().param_eta_glljqk[(g,l_star,l,j,q,k)]
+                    for l_star in b.parent_block().set_L[g]
+                    if l_star not in b.parent_block().set_L_exp[g]
+                    for j in b.parent_block().set_J[(g,l_star,l)]  # only directed arcs
+                ) == sum(b.var_trans_flows_s[s] for s in b.set_S)
+        b.constr_trans_flows = pyo.Constraint(rule=rule_constr_trans_flows)
+    
+        # *********************************************************************
+        # *********************************************************************
+        
+        # non-convex price functions
+        if not b.param_price_function_is_convex:
+            
+            # delta variables
+            b.var_active_segment_s = pyo.Var(b.set_S, within=pyo.Binary)
+            
+            # segments must be empty if the respective delta variable is zero
+            def rule_constr_empty_segment_if_delta_zero(b, s):
+                if len(b.set_S) == 1 or b.param_price_function_is_convex:
+                    # single segment, skip
+                    # convex, skip
+                    return pyo.Constraint.Skip
+                return (
+                    b.var_trans_flows_s[s] <= 
+                    b.param_v_max_s[s]*b.var_active_segment_s[s]
+                    )
+            b.constr_empty_segment_if_delta_zero = pyo.Constraint(
+                b.set_S, rule=rule_constr_empty_segment_if_delta_zero
+                )
+            
+            # if delta var is one, previous ones must be one too
+            # if delta var is zero, the next ones must also be zero
+            def rule_constr_delta_summing_logic(b, s):
+                if s == len(b.set_S)-1 or b.param_price_function_is_convex:
+                    # last segment, skip
+                    # convex, skip
+                    return pyo.Constraint.Skip
+                return (
+                    b.var_active_segment_s[s] >= 
+                    b.var_active_segment_s[s+1]
+                    )
+            b.constr_delta_summing_logic = pyo.Constraint(
+                b.set_S, rule=rule_constr_delta_summing_logic
+                )
+            
+            # if a segment is not completely used, the next ones must remain empty
+            def rule_constr_fill_up_segment_before_next(b, s):
+                if s == len(b.set_S)-1 or b.param_price_function_is_convex:
+                    # last segment, skip
+                    # convex, skip
+                    return pyo.Constraint.Skip
+                return (
+                    b.var_trans_flows_s[s] >= 
+                    b.var_active_segment_s[s+1]*
+                    b.param_v_max_s[s]
+                    )
+            b.constr_fill_up_segment_before_next = pyo.Constraint(
+                b.set_S, rule=rule_constr_fill_up_segment_before_next
+                )
+        
+        # *********************************************************************
+        # ********************************************************************* 
+        
+    model.block_prices = pyo.Block(model.set_GLQPK, rule=rule_block_prices)
+
+    # *************************************************************************
+    # *************************************************************************
+
+# *****************************************************************************
+# *****************************************************************************
+# *****************************************************************************
+# *****************************************************************************
+# *****************************************************************************
+# *****************************************************************************
+# *****************************************************************************
+# *****************************************************************************
+# *****************************************************************************
+# *****************************************************************************
+# *****************************************************************************
+# *****************************************************************************
+
+def price_lambda_no_block(
+    model: pyo.AbstractModel,
+    # convex_price_function: bool = True,
+    enable_default_values: bool = True,
+    enable_validation: bool = True,
+    enable_initialisation: bool = True
+    ):
+    
+    raise NotImplementedError
+    
+    # auxiliary set for pyomo
+    model.set_GLQPK = model.set_GL_exp_imp*model.set_QPK
+
+    # set of price segments
+    model.set_S = pyo.Set(model.set_GLQPK)
+
+    # set of GLQKS tuples
+    def init_set_GLQPKS(m):
+        return (
+            (g, l, q, p, k, s)
+            # for (g,l) in m.set_GL_exp_imp
+            # for (q,k) in m.set_QK
+            for (g, l, q, p, k) in m.set_S
+            for s in m.set_S[(g, l, q, p, k)]
+        )
+
+    model.set_GLQPKS = pyo.Set(
+        dimen=6, initialize=(init_set_GLQPKS if enable_initialisation else None)
+    )
+
+    # *************************************************************************
+    # *************************************************************************
+
+    # parameters
+
+    # resource prices
+
+    model.param_p_glqpks = pyo.Param(model.set_GLQPKS, within=pyo.NonNegativeReals)
+    
+    # price function convexity
+
+    model.param_price_function_is_convex = pyo.Param(
+        model.set_GLQPK, 
+        within=pyo.Boolean
+        )
+
+    # maximum resource volumes for each prices
+
+    model.param_v_max_glqpks = pyo.Param(
+        model.set_GLQPKS, 
+        within=pyo.NonNegativeReals,
+        # default=math.inf
+        )
+
+    # *************************************************************************
+    # *************************************************************************
+
+    # variables
+
+    # *************************************************************************
+    # *************************************************************************
+
+    # import and export flows
+    def bounds_var_trans_flows_glqpks(m, g, l, q, p, k, s):
+        # return (0, m.param_v_max_glqpks[(g, l, q, p, k, s)])
+        if (g, l, q, p, k, s) in m.param_v_max_glqpks:
+            # predefined finite capacity
+            return (0, m.param_v_max_glqpks[(g, l, q, p, k, s)])
+        else:
+            # infinite capacity
+            return (0, None)
+    model.var_trans_flows_glqpks = pyo.Var(
+        model.set_GLQPKS, within=pyo.NonNegativeReals, bounds=bounds_var_trans_flows_glqpks
+    )
+
+    # *************************************************************************
+    # *************************************************************************
+
+    # import flow costs and export flow revenues
+    def rule_constr_trans_monetary_flows(m, g, l, q, p, k):
+        if (g,l) in m.set_GL_imp:
+            return (
+                sum(
+                    m.var_trans_flows_glqpks[(g, l, q, p, k, s)]
+                    * m.param_p_glqpks[(g, l, q, p, k, s)]
+                    for s in m.set_S[(g, l, q, p, k)]
+                )
+                == m.var_ifc_glqpk[(g, l, q, p, k)]
+            )
+        else:
+            return (
+                sum(
+                    m.var_trans_flows_glqpks[(g, l, q, p, k, s)]
+                    * m.param_p_glqpks[(g, l, q, p, k, s)]
+                    for s in m.set_S[(g, l, q, p, k)]
+                )
+                == m.var_efr_glqpk[(g, l, q, p, k)]
+            )
+    model.constr_trans_monetary_flows = pyo.Constraint(
+        model.set_GLQPK, rule=rule_constr_trans_monetary_flows
+    )
+
+    # imported and exported flows
+    def rule_constr_trans_flows(m, g, l, q, p, k):
+        if (g,l) in m.set_GL_imp:
+            return sum(
+                m.var_v_glljqk[(g, l, l_star, j, q, k)]
+                for l_star in m.set_L[g]
+                if l_star not in m.set_L_imp[g]
+                for j in m.set_J[(g, l, l_star)]  # only directed arcs
+            ) == sum(m.var_trans_flows_glqpks[(g, l, q, p, k, s)] for s in m.set_S[(g, l, q, p, k)])
+        else:
+            return sum(
+                m.var_v_glljqk[(g, l_star, l, j, q, k)]
+                * m.param_eta_glljqk[(g, l_star, l, j, q, k)]
+                for l_star in m.set_L[g]
+                if l_star not in m.set_L_exp[g]
+                for j in m.set_J[(g, l_star, l)]  # only directed arcs
+            ) == sum(m.var_trans_flows_glqpks[(g, l, q, p, k, s)] for s in m.set_S[(g, l, q, p, k)])
+
+    model.constr_trans_flows = pyo.Constraint(
+        model.set_GLQPK, rule=rule_constr_trans_flows
+    )
+
+    # *************************************************************************
+    # *************************************************************************
+    
+    # non-convex price functions
+    # TODO: remove these variables from the model if they are not needed
+    # delta variables
+    model.var_active_segment_glqpks = pyo.Var(
+        model.set_GLQPKS, within=pyo.Binary
+    )
+    
+    # segments must be empty if the respective delta variable is zero
+    def rule_constr_empty_segment_if_delta_zero(m, g, l, q, p, k, s):
+        if len(m.set_S[(g,l,q,p,k)]) == 1 or m.param_price_function_is_convex[(g,l,q,p,k)]:
+            # single segment, skip
+            # convex, skip
+            return pyo.Constraint.Skip
+        return (
+            m.var_trans_flows_glqpks[(g,l,q,p,k,s)] <= 
+            m.param_v_max_glqpks[(g,l,q,p,k,s)]*
+            m.var_active_segment_glqpks[(g,l,q,p,k,s)]
+            )
+    model.constr_empty_segment_if_delta_zero = pyo.Constraint(
+        model.set_GLQPKS, rule=rule_constr_empty_segment_if_delta_zero
+        )
+    
+    # if delta var is one, previous ones must be one too
+    # if delta var is zero, the next ones must also be zero
+    def rule_constr_delta_summing_logic(m, g, l, q, p, k, s):
+        if s == len(m.set_S[(g,l,q,p,k)])-1 or m.param_price_function_is_convex[(g,l,q,p,k)]:
+            # last segment, skip
+            # convex, skip
+            return pyo.Constraint.Skip
+        return (
+            m.var_active_segment_glqpks[(g,l,q,p,k,s)] >= 
+            m.var_active_segment_glqpks[(g,l,q,p,k,s+1)]
+            )
+    model.constr_delta_summing_logic = pyo.Constraint(
+        model.set_GLQPKS, rule=rule_constr_delta_summing_logic
+        )
+    
+    # if a segment is not completely used, the next ones must remain empty
+    def rule_constr_fill_up_segment_before_next(m, g, l, q, p, k, s):
+        if s == len(m.set_S[(g,l,q,p,k)])-1 or m.param_price_function_is_convex[(g,l,q,p,k)]:
+            # last segment, skip
+            # convex, skip
+            return pyo.Constraint.Skip
+        return (
+            m.var_trans_flows_glqpks[(g,l,q,p,k,s)] >= 
+            m.var_active_segment_glqpks[(g,l,q,p,k,s+1)]*
+            m.param_v_max_glqpks[(g,l,q,p,k,s)]
+            )
+        # return (
+        #     m.var_if_glqpks[(g,l,q,p,k,s)]/m.param_v_max_glqpks[(g,l,q,p,k,s)] >= 
+        #     m.var_active_segment_glqpks[(g,l,q,p,k,s+1)]
+        #     )
+        # return (
+        #     m.param_v_max_glqpks[(g,l,q,p,k,s)]-m.var_if_glqpks[(g,l,q,p,k,s)] <= 
+        #     m.param_v_max_glqpks[(g,l,q,p,k,s)]*(1- m.var_active_segment_glqpks[(g,l,q,p,k,s+1)])
+        #     )
+    model.constr_fill_up_segment_before_next = pyo.Constraint(
+        model.set_GLQPKS, rule=rule_constr_fill_up_segment_before_next
+        )
+
+# *****************************************************************************
+# *****************************************************************************
\ No newline at end of file
diff --git a/src/topupopt/problems/esipp/blocks/other.py b/src/topupopt/problems/esipp/blocks/other.py
new file mode 100644
index 0000000000000000000000000000000000000000..ac4718040d3e7f90512278bba604a7f8c9586064
--- /dev/null
+++ b/src/topupopt/problems/esipp/blocks/other.py
@@ -0,0 +1,378 @@
+# imports
+
+import pyomo.environ as pyo
+
+# *****************************************************************************
+# *****************************************************************************
+
+# description: variation of the delta formulation
+
+# *****************************************************************************
+# *****************************************************************************
+
+def price_other_block(
+    model: pyo.AbstractModel,
+    enable_default_values: bool = True,
+    enable_validation: bool = True,
+    enable_initialisation: bool = True
+    ):
+    
+    # auxiliary set for pyomo
+    model.set_GLQPK = model.set_GL_exp_imp*model.set_QPK
+    
+    # create a block for a transshipment node during a given time interval
+    def rule_block_prices(b, g, l, q, p, k):
+        
+        # *********************************************************************
+        # *********************************************************************
+        
+        # sets
+        
+        # set of price segments
+        b.set_S = pyo.Set()
+        
+        # TODO: introduce a set of price segments for non-convex tariffs
+        # def init_set_S_nonconvex(b, s):
+        #     return (s for s in b.set_S if b.param_price_function_is_convex)
+        # b.set_S_nonconvex = pyo.Set(within=b.set_S, initialize=init_set_S_nonconvex)
+    
+        # *********************************************************************
+        # *********************************************************************
+        
+        # parameters
+
+        # resource prices
+        b.param_p_s = pyo.Param(b.set_S, within=pyo.NonNegativeReals)
+        
+        # price function convexity
+        b.param_price_function_is_convex = pyo.Param(within=pyo.Boolean)
+
+        # maximum resource volumes for each prices
+        b.param_v_max_s = pyo.Param(b.set_S, within=pyo.NonNegativeReals)
+    
+        # *********************************************************************
+        # *********************************************************************
+    
+        # variables
+    
+        # *********************************************************************
+        # *********************************************************************
+        
+        # TODO: consider replacing None in the bounds with inf in the parameter
+    
+        # import and export flows
+        def bounds_var_trans_flows_s(b, s):
+            if s in b.param_v_max_s:
+                # predefined finite capacity
+                return (0, b.param_v_max_s[s])
+            else:
+                # infinite capacity
+                return (0, None)
+        b.var_trans_flows_s = pyo.Var(
+            b.set_S, 
+            within=pyo.NonNegativeReals, 
+            bounds=bounds_var_trans_flows_s
+        )
+        
+        # *********************************************************************
+        # *********************************************************************
+    
+        # import flow costs and export flow revenues
+        def rule_constr_trans_monetary_flows(b):
+            if (g,l) in b.parent_block().set_GL_imp:
+                return (
+                    sum(b.var_trans_flows_s[s]*b.param_p_s[s]
+                        for s in b.set_S)
+                    == b.parent_block().var_ifc_glqpk[(g,l,q,p,k)]
+                )
+            else:
+                return (
+                    sum(b.var_trans_flows_s[s]*b.param_p_s[s]
+                        for s in b.set_S)
+                    == b.parent_block().var_efr_glqpk[(g,l,q,p,k)]
+                )
+        b.constr_trans_monetary_flows = pyo.Constraint(
+            rule=rule_constr_trans_monetary_flows
+        )
+        
+        # imported and exported flows
+        def rule_constr_trans_flows(b):
+            if (g,l) in b.parent_block().set_GL_imp:
+                return sum(
+                    b.parent_block().var_v_glljqk[(g,l,l_star,j,q,k)]
+                    for l_star in b.parent_block().set_L[g]
+                    if l_star not in b.parent_block().set_L_imp[g]
+                    for j in b.parent_block().set_J[(g,l,l_star)]  # only directed arcs
+                ) == sum(b.var_trans_flows_s[s] for s in b.set_S)
+            else:
+                return sum(
+                    b.parent_block().var_v_glljqk[(g,l_star,l,j,q,k)]
+                    * b.parent_block().param_eta_glljqk[(g,l_star,l,j,q,k)]
+                    for l_star in b.parent_block().set_L[g]
+                    if l_star not in b.parent_block().set_L_exp[g]
+                    for j in b.parent_block().set_J[(g,l_star,l)]  # only directed arcs
+                ) == sum(b.var_trans_flows_s[s] for s in b.set_S)
+        b.constr_trans_flows = pyo.Constraint(rule=rule_constr_trans_flows)
+    
+        # *********************************************************************
+        # *********************************************************************
+        
+        # non-convex price functions
+        if not b.param_price_function_is_convex:
+            
+            # delta variables
+            b.var_active_segment_s = pyo.Var(b.set_S, within=pyo.Binary)
+            
+            # segments must be empty if the respective delta variable is zero
+            def rule_constr_empty_segment_if_delta_zero(b, s):
+                if len(b.set_S) == 1 or b.param_price_function_is_convex:
+                    # single segment, skip
+                    # convex, skip
+                    return pyo.Constraint.Skip
+                return (
+                    b.var_trans_flows_s[s] <= 
+                    b.param_v_max_s[s]*b.var_active_segment_s[s]
+                    )
+            b.constr_empty_segment_if_delta_zero = pyo.Constraint(
+                b.set_S, rule=rule_constr_empty_segment_if_delta_zero
+                )
+            
+            # if delta var is one, previous ones must be one too
+            # if delta var is zero, the next ones must also be zero
+            def rule_constr_delta_summing_logic(b, s):
+                if s == len(b.set_S)-1 or b.param_price_function_is_convex:
+                    # last segment, skip
+                    # convex, skip
+                    return pyo.Constraint.Skip
+                return (
+                    b.var_active_segment_s[s] >= 
+                    b.var_active_segment_s[s+1]
+                    )
+            b.constr_delta_summing_logic = pyo.Constraint(
+                b.set_S, rule=rule_constr_delta_summing_logic
+                )
+            
+            # if a segment is not completely used, the next ones must remain empty
+            def rule_constr_fill_up_segment_before_next(b, s):
+                if s == len(b.set_S)-1 or b.param_price_function_is_convex:
+                    # last segment, skip
+                    # convex, skip
+                    return pyo.Constraint.Skip
+                return (
+                    b.var_trans_flows_s[s] >= 
+                    b.var_active_segment_s[s+1]*
+                    b.param_v_max_s[s]
+                    )
+            b.constr_fill_up_segment_before_next = pyo.Constraint(
+                b.set_S, rule=rule_constr_fill_up_segment_before_next
+                )
+        
+        # *********************************************************************
+        # ********************************************************************* 
+        
+    model.block_prices = pyo.Block(model.set_GLQPK, rule=rule_block_prices)
+
+    # *************************************************************************
+    # *************************************************************************
+    
+
+# *****************************************************************************
+# *****************************************************************************
+# *****************************************************************************
+# *****************************************************************************
+# *****************************************************************************
+# *****************************************************************************
+# *****************************************************************************
+# *****************************************************************************
+# *****************************************************************************
+# *****************************************************************************
+# *****************************************************************************
+# *****************************************************************************
+
+def price_other_no_block(
+    model: pyo.AbstractModel,
+    # convex_price_function: bool = True,
+    enable_default_values: bool = True,
+    enable_validation: bool = True,
+    enable_initialisation: bool = True
+    ):
+    
+    # auxiliary set for pyomo
+    model.set_GLQPK = model.set_GL_exp_imp*model.set_QPK
+
+    # set of price segments
+    model.set_S = pyo.Set(model.set_GLQPK)
+
+    # set of GLQKS tuples
+    def init_set_GLQPKS(m):
+        return (
+            (g, l, q, p, k, s)
+            # for (g,l) in m.set_GL_exp_imp
+            # for (q,k) in m.set_QK
+            for (g, l, q, p, k) in m.set_S
+            for s in m.set_S[(g, l, q, p, k)]
+        )
+
+    model.set_GLQPKS = pyo.Set(
+        dimen=6, initialize=(init_set_GLQPKS if enable_initialisation else None)
+    )
+
+    # *************************************************************************
+    # *************************************************************************
+
+    # parameters
+
+    # resource prices
+
+    model.param_p_glqpks = pyo.Param(model.set_GLQPKS, within=pyo.NonNegativeReals)
+    
+    # price function convexity
+
+    model.param_price_function_is_convex = pyo.Param(
+        model.set_GLQPK, 
+        within=pyo.Boolean
+        )
+
+    # maximum resource volumes for each prices
+
+    model.param_v_max_glqpks = pyo.Param(
+        model.set_GLQPKS, 
+        within=pyo.NonNegativeReals,
+        # default=math.inf
+        )
+
+    # *************************************************************************
+    # *************************************************************************
+
+    # variables
+
+    # *************************************************************************
+    # *************************************************************************
+
+    # import and export flows
+    def bounds_var_trans_flows_glqpks(m, g, l, q, p, k, s):
+        # return (0, m.param_v_max_glqpks[(g, l, q, p, k, s)])
+        if (g, l, q, p, k, s) in m.param_v_max_glqpks:
+            # predefined finite capacity
+            return (0, m.param_v_max_glqpks[(g, l, q, p, k, s)])
+        else:
+            # infinite capacity
+            return (0, None)
+    model.var_trans_flows_glqpks = pyo.Var(
+        model.set_GLQPKS, within=pyo.NonNegativeReals, bounds=bounds_var_trans_flows_glqpks
+    )
+
+    # *************************************************************************
+    # *************************************************************************
+
+    # import flow costs and export flow revenues
+    def rule_constr_trans_monetary_flows(m, g, l, q, p, k):
+        if (g,l) in m.set_GL_imp:
+            return (
+                sum(
+                    m.var_trans_flows_glqpks[(g, l, q, p, k, s)]
+                    * m.param_p_glqpks[(g, l, q, p, k, s)]
+                    for s in m.set_S[(g, l, q, p, k)]
+                )
+                == m.var_ifc_glqpk[(g, l, q, p, k)]
+            )
+        else:
+            return (
+                sum(
+                    m.var_trans_flows_glqpks[(g, l, q, p, k, s)]
+                    * m.param_p_glqpks[(g, l, q, p, k, s)]
+                    for s in m.set_S[(g, l, q, p, k)]
+                )
+                == m.var_efr_glqpk[(g, l, q, p, k)]
+            )
+    model.constr_trans_monetary_flows = pyo.Constraint(
+        model.set_GLQPK, rule=rule_constr_trans_monetary_flows
+    )
+
+    # imported and exported flows
+    def rule_constr_trans_flows(m, g, l, q, p, k):
+        if (g,l) in m.set_GL_imp:
+            return sum(
+                m.var_v_glljqk[(g, l, l_star, j, q, k)]
+                for l_star in m.set_L[g]
+                if l_star not in m.set_L_imp[g]
+                for j in m.set_J[(g, l, l_star)]  # only directed arcs
+            ) == sum(m.var_trans_flows_glqpks[(g, l, q, p, k, s)] for s in m.set_S[(g, l, q, p, k)])
+        else:
+            return sum(
+                m.var_v_glljqk[(g, l_star, l, j, q, k)]
+                * m.param_eta_glljqk[(g, l_star, l, j, q, k)]
+                for l_star in m.set_L[g]
+                if l_star not in m.set_L_exp[g]
+                for j in m.set_J[(g, l_star, l)]  # only directed arcs
+            ) == sum(m.var_trans_flows_glqpks[(g, l, q, p, k, s)] for s in m.set_S[(g, l, q, p, k)])
+
+    model.constr_trans_flows = pyo.Constraint(
+        model.set_GLQPK, rule=rule_constr_trans_flows
+    )
+
+    # *************************************************************************
+    # *************************************************************************
+    
+    # non-convex price functions
+    # TODO: remove these variables from the model if they are not needed
+    # delta variables
+    model.var_active_segment_glqpks = pyo.Var(
+        model.set_GLQPKS, within=pyo.Binary
+    )
+    
+    # segments must be empty if the respective delta variable is zero
+    def rule_constr_empty_segment_if_delta_zero(m, g, l, q, p, k, s):
+        if len(m.set_S[(g,l,q,p,k)]) == 1 or m.param_price_function_is_convex[(g,l,q,p,k)]:
+            # single segment, skip
+            # convex, skip
+            return pyo.Constraint.Skip
+        return (
+            m.var_trans_flows_glqpks[(g,l,q,p,k,s)] <= 
+            m.param_v_max_glqpks[(g,l,q,p,k,s)]*
+            m.var_active_segment_glqpks[(g,l,q,p,k,s)]
+            )
+    model.constr_empty_segment_if_delta_zero = pyo.Constraint(
+        model.set_GLQPKS, rule=rule_constr_empty_segment_if_delta_zero
+        )
+    
+    # if delta var is one, previous ones must be one too
+    # if delta var is zero, the next ones must also be zero
+    def rule_constr_delta_summing_logic(m, g, l, q, p, k, s):
+        if s == len(m.set_S[(g,l,q,p,k)])-1 or m.param_price_function_is_convex[(g,l,q,p,k)]:
+            # last segment, skip
+            # convex, skip
+            return pyo.Constraint.Skip
+        return (
+            m.var_active_segment_glqpks[(g,l,q,p,k,s)] >= 
+            m.var_active_segment_glqpks[(g,l,q,p,k,s+1)]
+            )
+    model.constr_delta_summing_logic = pyo.Constraint(
+        model.set_GLQPKS, rule=rule_constr_delta_summing_logic
+        )
+    
+    # if a segment is not completely used, the next ones must remain empty
+    def rule_constr_fill_up_segment_before_next(m, g, l, q, p, k, s):
+        if s == len(m.set_S[(g,l,q,p,k)])-1 or m.param_price_function_is_convex[(g,l,q,p,k)]:
+            # last segment, skip
+            # convex, skip
+            return pyo.Constraint.Skip
+        return (
+            m.var_trans_flows_glqpks[(g,l,q,p,k,s)] >= 
+            m.var_active_segment_glqpks[(g,l,q,p,k,s+1)]*
+            m.param_v_max_glqpks[(g,l,q,p,k,s)]
+            )
+        # return (
+        #     m.var_if_glqpks[(g,l,q,p,k,s)]/m.param_v_max_glqpks[(g,l,q,p,k,s)] >= 
+        #     m.var_active_segment_glqpks[(g,l,q,p,k,s+1)]
+        #     )
+        # return (
+        #     m.param_v_max_glqpks[(g,l,q,p,k,s)]-m.var_if_glqpks[(g,l,q,p,k,s)] <= 
+        #     m.param_v_max_glqpks[(g,l,q,p,k,s)]*(1- m.var_active_segment_glqpks[(g,l,q,p,k,s+1)])
+        #     )
+    model.constr_fill_up_segment_before_next = pyo.Constraint(
+        model.set_GLQPKS, rule=rule_constr_fill_up_segment_before_next
+        )
+
+# *****************************************************************************
+# *****************************************************************************
\ No newline at end of file
diff --git a/src/topupopt/problems/esipp/blocks/prices.py b/src/topupopt/problems/esipp/blocks/prices.py
index 097e9356cbb2541b84513690ea91fd50475fef8f..4c38b3fe092c8b0d23b56cb0030acc86e52ddd87 100644
--- a/src/topupopt/problems/esipp/blocks/prices.py
+++ b/src/topupopt/problems/esipp/blocks/prices.py
@@ -1,14 +1,25 @@
 # imports
 
 import pyomo.environ as pyo
-import math
+from .other import price_other_block, price_other_no_block
+from .delta import price_delta_block, price_delta_no_block
 
 # *****************************************************************************
 # *****************************************************************************
 
+NODE_PRICE_LAMBDA = 'lambda'
+NODE_PRICE_DELTA = 'delta'
+NODE_PRICE_OTHER = None
+NODE_PRICES = (
+    NODE_PRICE_LAMBDA,
+    NODE_PRICE_DELTA,
+    NODE_PRICE_OTHER
+    )
+
 def add_price_functions(
     model: pyo.AbstractModel,
-    use_blocks: bool = False,
+    use_blocks: bool = True,
+    node_price_model = NODE_PRICE_OTHER,
     **kwargs
 ):  
     
@@ -16,399 +27,24 @@ def add_price_functions(
     # *************************************************************************
     
     if use_blocks:
-        # blocks
-        price_other_block(model, **kwargs)
-    else:
-        # no blocks
-        price_other(model, **kwargs)
-        
-    # *************************************************************************
-    # *************************************************************************
-
-# *****************************************************************************
-# *****************************************************************************
-
-def price_other_block(
-    model: pyo.AbstractModel,
-    # convex_price_function: bool = True,
-    enable_default_values: bool = True,
-    enable_validation: bool = True,
-    enable_initialisation: bool = True
-    ):
-    
-    # auxiliary set for pyomo
-    model.set_GLQPK = model.set_GL_exp_imp*model.set_QPK
-    
-    # create a block for a transshipment node during a given time interval
-    def rule_block_prices(b, g, l, q, p, k):
-        
-        # *********************************************************************
-        # *********************************************************************
-        
-        # sets
-        
-        # set of price segments
-        b.set_S = pyo.Set()
-        
-        # TODO: introduce a set of price segments for non-convex tariffs
-        # def init_set_S_nonconvex(b, s):
-        #     return (s for s in b.set_S if b.param_price_function_is_convex)
-        # b.set_S_nonconvex = pyo.Set(within=b.set_S, initialize=init_set_S_nonconvex)
-    
-        # *********************************************************************
-        # *********************************************************************
-        
-        # parameters
-
-        # resource prices
-        b.param_p_s = pyo.Param(b.set_S, within=pyo.NonNegativeReals)
-        
-        # price function convexity
-        b.param_price_function_is_convex = pyo.Param(within=pyo.Boolean)
-
-        # maximum resource volumes for each prices
-        b.param_v_max_s = pyo.Param(b.set_S, within=pyo.NonNegativeReals)
-    
-        # *********************************************************************
-        # *********************************************************************
-    
-        # variables
-    
-        # *********************************************************************
-        # *********************************************************************
-        
-        # TODO: consider replacing None in the bounds with inf in the parameter
-    
-        # import and export flows
-        def bounds_var_trans_flows_s(b, s):
-            if s in b.param_v_max_s:
-                # predefined finite capacity
-                return (0, b.param_v_max_s[s])
-            else:
-                # infinite capacity
-                return (0, None)
-        b.var_trans_flows_s = pyo.Var(
-            b.set_S, 
-            within=pyo.NonNegativeReals, 
-            bounds=bounds_var_trans_flows_s
-        )
-        
-        # *********************************************************************
-        # *********************************************************************
-    
-        # import flow costs and export flow revenues
-        def rule_constr_trans_monetary_flows(b):
-            if (g,l) in b.parent_block().set_GL_imp:
-                return (
-                    sum(b.var_trans_flows_s[s]*b.param_p_s[s]
-                        for s in b.set_S)
-                    == b.parent_block().var_ifc_glqpk[(g,l,q,p,k)]
-                )
-            else:
-                return (
-                    sum(b.var_trans_flows_s[s]*b.param_p_s[s]
-                        for s in b.set_S)
-                    == b.parent_block().var_efr_glqpk[(g,l,q,p,k)]
-                )
-        b.constr_trans_monetary_flows = pyo.Constraint(
-            rule=rule_constr_trans_monetary_flows
-        )
-        
-        # imported and exported flows
-        def rule_constr_trans_flows(b):
-            if (g,l) in b.parent_block().set_GL_imp:
-                return sum(
-                    b.parent_block().var_v_glljqk[(g,l,l_star,j,q,k)]
-                    for l_star in b.parent_block().set_L[g]
-                    if l_star not in b.parent_block().set_L_imp[g]
-                    for j in b.parent_block().set_J[(g,l,l_star)]  # only directed arcs
-                ) == sum(b.var_trans_flows_s[s] for s in b.set_S)
-            else:
-                return sum(
-                    b.parent_block().var_v_glljqk[(g,l_star,l,j,q,k)]
-                    * b.parent_block().param_eta_glljqk[(g,l_star,l,j,q,k)]
-                    for l_star in b.parent_block().set_L[g]
-                    if l_star not in b.parent_block().set_L_exp[g]
-                    for j in b.parent_block().set_J[(g,l_star,l)]  # only directed arcs
-                ) == sum(b.var_trans_flows_s[s] for s in b.set_S)
-        b.constr_trans_flows = pyo.Constraint(rule=rule_constr_trans_flows)
-    
-        # *********************************************************************
-        # *********************************************************************
-        
-        # non-convex price functions
-        if not b.param_price_function_is_convex:
-            
-            # delta variables
-            b.var_active_segment_s = pyo.Var(b.set_S, within=pyo.Binary)
-            
-            # segments must be empty if the respective delta variable is zero
-            def rule_constr_empty_segment_if_delta_zero(b, s):
-                if len(b.set_S) == 1 or b.param_price_function_is_convex:
-                    # single segment, skip
-                    # convex, skip
-                    return pyo.Constraint.Skip
-                return (
-                    b.var_trans_flows_s[s] <= 
-                    b.param_v_max_s[s]*b.var_active_segment_s[s]
-                    )
-            b.constr_empty_segment_if_delta_zero = pyo.Constraint(
-                b.set_S, rule=rule_constr_empty_segment_if_delta_zero
-                )
-            
-            # if delta var is one, previous ones must be one too
-            # if delta var is zero, the next ones must also be zero
-            def rule_constr_delta_summing_logic(b, s):
-                if s == len(b.set_S)-1 or b.param_price_function_is_convex:
-                    # last segment, skip
-                    # convex, skip
-                    return pyo.Constraint.Skip
-                return (
-                    b.var_active_segment_s[s] >= 
-                    b.var_active_segment_s[s+1]
-                    )
-            b.constr_delta_summing_logic = pyo.Constraint(
-                b.set_S, rule=rule_constr_delta_summing_logic
-                )
-            
-            # if a segment is not completely used, the next ones must remain empty
-            def rule_constr_fill_up_segment_before_next(b, s):
-                if s == len(b.set_S)-1 or b.param_price_function_is_convex:
-                    # last segment, skip
-                    # convex, skip
-                    return pyo.Constraint.Skip
-                return (
-                    b.var_trans_flows_s[s] >= 
-                    b.var_active_segment_s[s+1]*
-                    b.param_v_max_s[s]
-                    )
-            b.constr_fill_up_segment_before_next = pyo.Constraint(
-                b.set_S, rule=rule_constr_fill_up_segment_before_next
-                )
-        
-        # *********************************************************************
-        # ********************************************************************* 
-        
-    model.block_prices = pyo.Block(model.set_GLQPK, rule=rule_block_prices)
-
-    # *************************************************************************
-    # *************************************************************************
-
-# *****************************************************************************
-# *****************************************************************************
-# *****************************************************************************
-# *****************************************************************************
-# *****************************************************************************
-# *****************************************************************************
-# *****************************************************************************
-# *****************************************************************************
-# *****************************************************************************
-# *****************************************************************************
-# *****************************************************************************
-# *****************************************************************************
-
-def price_other(
-    model: pyo.AbstractModel,
-    # convex_price_function: bool = True,
-    enable_default_values: bool = True,
-    enable_validation: bool = True,
-    enable_initialisation: bool = True
-    ):
-    
-    # auxiliary set for pyomo
-    model.set_GLQPK = model.set_GL_exp_imp*model.set_QPK
-
-    # set of price segments
-    model.set_S = pyo.Set(model.set_GLQPK)
-
-    # set of GLQKS tuples
-    def init_set_GLQPKS(m):
-        return (
-            (g, l, q, p, k, s)
-            # for (g,l) in m.set_GL_exp_imp
-            # for (q,k) in m.set_QK
-            for (g, l, q, p, k) in m.set_S
-            for s in m.set_S[(g, l, q, p, k)]
-        )
-
-    model.set_GLQPKS = pyo.Set(
-        dimen=6, initialize=(init_set_GLQPKS if enable_initialisation else None)
-    )
-
-    # *************************************************************************
-    # *************************************************************************
-
-    # parameters
-
-    # resource prices
-
-    model.param_p_glqpks = pyo.Param(model.set_GLQPKS, within=pyo.NonNegativeReals)
-    
-    # price function convexity
-
-    model.param_price_function_is_convex = pyo.Param(
-        model.set_GLQPK, 
-        within=pyo.Boolean
-        )
-
-    # maximum resource volumes for each prices
-
-    model.param_v_max_glqpks = pyo.Param(
-        model.set_GLQPKS, 
-        within=pyo.NonNegativeReals,
-        # default=math.inf
-        )
-
-    # *************************************************************************
-    # *************************************************************************
-
-    # variables
-
-    # *************************************************************************
-    # *************************************************************************
-
-    # import and export flows
-    def bounds_var_trans_flows_glqpks(m, g, l, q, p, k, s):
-        # return (0, m.param_v_max_glqpks[(g, l, q, p, k, s)])
-        if (g, l, q, p, k, s) in m.param_v_max_glqpks:
-            # predefined finite capacity
-            return (0, m.param_v_max_glqpks[(g, l, q, p, k, s)])
-        else:
-            # infinite capacity
-            return (0, None)
-    model.var_trans_flows_glqpks = pyo.Var(
-        model.set_GLQPKS, within=pyo.NonNegativeReals, bounds=bounds_var_trans_flows_glqpks
-    )
-
-    # *************************************************************************
-    # *************************************************************************
-
-    # import flow costs and export flow revenues
-    def rule_constr_trans_monetary_flows(m, g, l, q, p, k):
-        if (g,l) in m.set_GL_imp:
-            return (
-                sum(
-                    m.var_trans_flows_glqpks[(g, l, q, p, k, s)]
-                    * m.param_p_glqpks[(g, l, q, p, k, s)]
-                    for s in m.set_S[(g, l, q, p, k)]
-                )
-                == m.var_ifc_glqpk[(g, l, q, p, k)]
-            )
+        # with blocks
+        if node_price_model == NODE_PRICE_LAMBDA:
+            raise NotImplementedError 
+        elif node_price_model == NODE_PRICE_DELTA:
+            price_delta_block(model, **kwargs)
         else:
-            return (
-                sum(
-                    m.var_trans_flows_glqpks[(g, l, q, p, k, s)]
-                    * m.param_p_glqpks[(g, l, q, p, k, s)]
-                    for s in m.set_S[(g, l, q, p, k)]
-                )
-                == m.var_efr_glqpk[(g, l, q, p, k)]
-            )
-    model.constr_trans_monetary_flows = pyo.Constraint(
-        model.set_GLQPK, rule=rule_constr_trans_monetary_flows
-    )
-
-    # imported and exported flows
-    def rule_constr_trans_flows(m, g, l, q, p, k):
-        if (g,l) in m.set_GL_imp:
-            return sum(
-                m.var_v_glljqk[(g, l, l_star, j, q, k)]
-                for l_star in m.set_L[g]
-                if l_star not in m.set_L_imp[g]
-                for j in m.set_J[(g, l, l_star)]  # only directed arcs
-            ) == sum(m.var_trans_flows_glqpks[(g, l, q, p, k, s)] for s in m.set_S[(g, l, q, p, k)])
+            price_other_block(model, **kwargs)
+    else:
+        # without blocks
+        if node_price_model == NODE_PRICE_LAMBDA:
+            raise NotImplementedError 
+        elif node_price_model == NODE_PRICE_DELTA:
+            raise NotImplementedError 
         else:
-            return sum(
-                m.var_v_glljqk[(g, l_star, l, j, q, k)]
-                * m.param_eta_glljqk[(g, l_star, l, j, q, k)]
-                for l_star in m.set_L[g]
-                if l_star not in m.set_L_exp[g]
-                for j in m.set_J[(g, l_star, l)]  # only directed arcs
-            ) == sum(m.var_trans_flows_glqpks[(g, l, q, p, k, s)] for s in m.set_S[(g, l, q, p, k)])
-
-    model.constr_trans_flows = pyo.Constraint(
-        model.set_GLQPK, rule=rule_constr_trans_flows
-    )
-
+            price_other_no_block(model, **kwargs)
+        
     # *************************************************************************
     # *************************************************************************
-    
-    # non-convex price functions
-    # TODO: remove these variables from the model if they are not needed
-    # delta variables
-    model.var_active_segment_glqpks = pyo.Var(
-        model.set_GLQPKS, within=pyo.Binary
-    )
-    
-    # segments must be empty if the respective delta variable is zero
-    def rule_constr_empty_segment_if_delta_zero(m, g, l, q, p, k, s):
-        if len(m.set_S[(g,l,q,p,k)]) == 1 or m.param_price_function_is_convex[(g,l,q,p,k)]:
-            # single segment, skip
-            # convex, skip
-            return pyo.Constraint.Skip
-        return (
-            m.var_trans_flows_glqpks[(g,l,q,p,k,s)] <= 
-            m.param_v_max_glqpks[(g,l,q,p,k,s)]*
-            m.var_active_segment_glqpks[(g,l,q,p,k,s)]
-            )
-    model.constr_empty_segment_if_delta_zero = pyo.Constraint(
-        model.set_GLQPKS, rule=rule_constr_empty_segment_if_delta_zero
-        )
-    
-    # if delta var is one, previous ones must be one too
-    # if delta var is zero, the next ones must also be zero
-    def rule_constr_delta_summing_logic(m, g, l, q, p, k, s):
-        if s == len(m.set_S[(g,l,q,p,k)])-1 or m.param_price_function_is_convex[(g,l,q,p,k)]:
-            # last segment, skip
-            # convex, skip
-            return pyo.Constraint.Skip
-        return (
-            m.var_active_segment_glqpks[(g,l,q,p,k,s)] >= 
-            m.var_active_segment_glqpks[(g,l,q,p,k,s+1)]
-            )
-    model.constr_delta_summing_logic = pyo.Constraint(
-        model.set_GLQPKS, rule=rule_constr_delta_summing_logic
-        )
-    
-    # if a segment is not completely used, the next ones must remain empty
-    def rule_constr_fill_up_segment_before_next(m, g, l, q, p, k, s):
-        if s == len(m.set_S[(g,l,q,p,k)])-1 or m.param_price_function_is_convex[(g,l,q,p,k)]:
-            # last segment, skip
-            # convex, skip
-            return pyo.Constraint.Skip
-        return (
-            m.var_trans_flows_glqpks[(g,l,q,p,k,s)] >= 
-            m.var_active_segment_glqpks[(g,l,q,p,k,s+1)]*
-            m.param_v_max_glqpks[(g,l,q,p,k,s)]
-            )
-        # return (
-        #     m.var_if_glqpks[(g,l,q,p,k,s)]/m.param_v_max_glqpks[(g,l,q,p,k,s)] >= 
-        #     m.var_active_segment_glqpks[(g,l,q,p,k,s+1)]
-        #     )
-        # return (
-        #     m.param_v_max_glqpks[(g,l,q,p,k,s)]-m.var_if_glqpks[(g,l,q,p,k,s)] <= 
-        #     m.param_v_max_glqpks[(g,l,q,p,k,s)]*(1- m.var_active_segment_glqpks[(g,l,q,p,k,s+1)])
-        #     )
-    model.constr_fill_up_segment_before_next = pyo.Constraint(
-        model.set_GLQPKS, rule=rule_constr_fill_up_segment_before_next
-        )
-
-# *****************************************************************************
-# *****************************************************************************
-
-# TODO: implement the lambda model
-
-def price_block_lambda(model: pyo.AbstractModel, **kwargs):
-    
-    raise NotImplementedError
-
-# *****************************************************************************
-# *****************************************************************************
-
-# TODO: implement the delta model
-
-def price_block_delta(model: pyo.AbstractModel, **kwargs):
-    
-    raise NotImplementedError
 
 # *****************************************************************************
 # *****************************************************************************
\ No newline at end of file
diff --git a/src/topupopt/problems/esipp/model.py b/src/topupopt/problems/esipp/model.py
index 508ed3e9ddf2c71b4d1f277de99e8889db205e4c..68bb0b3f2a80207f5f06945841176f0ffd536ee5 100644
--- a/src/topupopt/problems/esipp/model.py
+++ b/src/topupopt/problems/esipp/model.py
@@ -11,10 +11,11 @@ from .blocks.prices import add_price_functions
 
 def create_model(
     name: str,
+    use_prices_block: bool,
+    node_price_model: str,
     enable_default_values: bool = True,
     enable_validation: bool = True,
-    enable_initialisation: bool = True,
-    use_prices_block: bool = False
+    enable_initialisation: bool = True
 ):
     # TODO: test initialisation
 
@@ -2131,7 +2132,11 @@ def create_model(
     )
     
     # prices
-    add_price_functions(model, use_blocks=use_prices_block)
+    add_price_functions(
+        model,
+        use_blocks=use_prices_block, 
+        node_price_model=node_price_model
+        )
 
     # *************************************************************************
     # *************************************************************************
diff --git a/src/topupopt/problems/esipp/problem.py b/src/topupopt/problems/esipp/problem.py
index 2770f9c2a39c480895955e9955115cecdca903a1..4d89b122ee5940d4a3c7a56aa7e2a35cbb938221 100644
--- a/src/topupopt/problems/esipp/problem.py
+++ b/src/topupopt/problems/esipp/problem.py
@@ -16,6 +16,7 @@ from .system import EnergySystem
 from .resource import ResourcePrice
 from .time import EconomicTimeFrame
 from .time import entries_are_invariant
+from .blocks.prices import NODE_PRICE_OTHER, NODE_PRICES
 
 # *****************************************************************************
 # *****************************************************************************
@@ -63,15 +64,6 @@ class InfrastructurePlanningProblem(EnergySystem):
         STATIC_LOSS_MODE_US,
         STATIC_LOSS_MODE_DS,
     )
-    
-    NODE_PRICE_LAMBDA = 1
-    NODE_PRICE_DELTA = 2
-    NODE_PRICE_OTHER = 3
-    NODE_PRICES = (
-        NODE_PRICE_LAMBDA,
-        NODE_PRICE_DELTA,
-        NODE_PRICE_OTHER
-        )
 
     # *************************************************************************
     # *************************************************************************
@@ -89,7 +81,7 @@ class InfrastructurePlanningProblem(EnergySystem):
         converters: dict = None,
         prepare_model: bool = True,
         validate_inputs: bool = True,
-        node_price_model = NODE_PRICE_DELTA,
+        node_price_model = NODE_PRICE_OTHER,
         use_prices_block: bool = True #False
     ):  # TODO: switch to False when everything is more mature
         # *********************************************************************
@@ -129,6 +121,7 @@ class InfrastructurePlanningProblem(EnergySystem):
         
         # options
         self.use_prices_block = use_prices_block
+        self.node_price_model = node_price_model
         
         # *********************************************************************
         # *********************************************************************
@@ -1668,7 +1661,11 @@ class InfrastructurePlanningProblem(EnergySystem):
     def prepare(self, name = None):
         """Sets up the problem model with which instances can be built."""
         # create pyomo model (AbstractModel)
-        self.model = create_model(name, use_prices_block=self.use_prices_block)
+        self.model = create_model(
+            name, 
+            use_prices_block=self.use_prices_block,
+            node_price_model=self.node_price_model
+            )
 
     # *************************************************************************
     # *************************************************************************
diff --git a/src/topupopt/problems/esipp/utils.py b/src/topupopt/problems/esipp/utils.py
index 08f6f612f55dae9f730c097dbfb027f6dd812fa1..a61b35a9ac94ad5ae7228308ab8baeb8703d7542 100644
--- a/src/topupopt/problems/esipp/utils.py
+++ b/src/topupopt/problems/esipp/utils.py
@@ -13,6 +13,7 @@ from matplotlib import pyplot as plt
 # local, internal
 from .problem import InfrastructurePlanningProblem
 from .network import Network
+from .blocks.prices import NODE_PRICE_DELTA, NODE_PRICE_LAMBDA #, NODE_PRICE_OTHER
 
 # *****************************************************************************
 # *****************************************************************************
@@ -97,34 +98,58 @@ def statistics(ipp: InfrastructurePlanningProblem,
             )
     if prices_in_block:
         # imports
-        imports_qpk = {
-            qpk: pyo.value( 
-                sum(
-                    ipp.instance.block_prices[(g,l_imp,*qpk)].var_trans_flows_s[s]
-                    for g, l_imp in import_node_keys
-                    # for g in ipp.networks
-                    # for l_imp in ipp.networks[g].import_nodes
-                    for s in ipp.instance.block_prices[(g,l_imp,*qpk)].set_S
+        if ipp.node_price_model == NODE_PRICE_DELTA:
+            # imports
+            imports_qpk = {
+                qpk: pyo.value( 
+                    sum(
+                        ipp.instance.block_prices[(g,l_imp,*qpk)].var_trans_flows
+                        for g, l_imp in import_node_keys
+                        )
+                    *ipp.instance.param_c_time_qpk[qpk]
                     )
-                *ipp.instance.param_c_time_qpk[qpk]
-                )
-            for qpk in ipp.time_frame.qpk()
-            }
-        
-        # exports
-        exports_qpk = {
-            qpk: pyo.value(
-                sum(
-                    ipp.instance.block_prices[(g,l_exp,*qpk)].var_trans_flows_s[s]
-                    for g, l_exp in export_node_keys
-                    # for g in ipp.networks
-                    # for l_exp in ipp.networks[g].export_nodes
-                    for s in ipp.instance.block_prices[(g,l_exp,*qpk)].set_S
+                for qpk in ipp.time_frame.qpk()
+                }
+            # exports
+            exports_qpk = {
+                qpk: pyo.value(
+                    sum(
+                        ipp.instance.block_prices[(g,l_exp,*qpk)].var_trans_flows
+                        for g, l_exp in export_node_keys
+                        )
+                    *ipp.instance.param_c_time_qpk[qpk]
                     )
-                *ipp.instance.param_c_time_qpk[qpk]
-                )
-            for qpk in ipp.time_frame.qpk()
-            }
+                for qpk in ipp.time_frame.qpk()
+                }
+        else:
+            # imports
+            imports_qpk = {
+                qpk: pyo.value( 
+                    sum(
+                        ipp.instance.block_prices[(g,l_imp,*qpk)].var_trans_flows_s[s]
+                        for g, l_imp in import_node_keys
+                        # for g in ipp.networks
+                        # for l_imp in ipp.networks[g].import_nodes
+                        for s in ipp.instance.block_prices[(g,l_imp,*qpk)].set_S
+                        )
+                    *ipp.instance.param_c_time_qpk[qpk]
+                    )
+                for qpk in ipp.time_frame.qpk()
+                }
+            # exports
+            exports_qpk = {
+                qpk: pyo.value(
+                    sum(
+                        ipp.instance.block_prices[(g,l_exp,*qpk)].var_trans_flows_s[s]
+                        for g, l_exp in export_node_keys
+                        # for g in ipp.networks
+                        # for l_exp in ipp.networks[g].export_nodes
+                        for s in ipp.instance.block_prices[(g,l_exp,*qpk)].set_S
+                        )
+                    *ipp.instance.param_c_time_qpk[qpk]
+                    )
+                for qpk in ipp.time_frame.qpk()
+                }
     else:
         # not in a block
         # imports
diff --git a/tests/test_esipp.py b/tests/test_esipp.py
index d134a51c71f78474ad0dcb1ce02aa37f775fbf06..570c6edb5b1ecc2efbf51ea64aba2037b00f8f97 100644
--- a/tests/test_esipp.py
+++ b/tests/test_esipp.py
@@ -6,6 +6,7 @@
 from src.topupopt.problems.esipp.problem import InfrastructurePlanningProblem
 from src.topupopt.problems.esipp.network import Network
 from src.topupopt.problems.esipp.time import EconomicTimeFrame
+from src.topupopt.problems.esipp.blocks.prices import NODE_PRICE_OTHER
 
 # *****************************************************************************
 # *****************************************************************************
@@ -46,7 +47,8 @@ def build_solve_ipp(
     # discount_rates: dict = None,
     assessment_weights: dict = None,
     simplify_problem: bool = False,
-    use_prices_block: bool = False
+    use_prices_block: bool = False,
+    node_price_model: int = NODE_PRICE_OTHER
 ):
     
     if type(assessment_weights) != dict:
@@ -81,7 +83,8 @@ def build_solve_ipp(
         time_weights=time_weights,
         normalised_time_interval_duration=normalised_time_interval_duration,
         assessment_weights=assessment_weights,
-        use_prices_block=use_prices_block
+        use_prices_block=use_prices_block,
+        node_price_model=node_price_model
     )
 
     # add networks and systems
@@ -222,7 +225,7 @@ def build_solve_ipp(
     # ipp.instantiate(place_fixed_losses_upstream_if_possible=False)
 
     ipp.instantiate(initialise_ancillary_sets=init_aux_sets)
-    # ipp.instance.pprint()
+    
     # optimise
     ipp.optimise(
         solver_name=solver,
diff --git a/tests/test_esipp_prices.py b/tests/test_esipp_prices.py
index 704e488a802984df87b74d60fd7e7865dd2c424f..8cd73ddcfb2e7279fb250e8b51975c92ba3dfd26 100644
--- a/tests/test_esipp_prices.py
+++ b/tests/test_esipp_prices.py
@@ -20,6 +20,7 @@ from src.topupopt.problems.esipp.utils import statistics
 from src.topupopt.problems.esipp.time import EconomicTimeFrame
 # from src.topupopt.problems.esipp.converter import Converter
 from test_esipp import build_solve_ipp
+from src.topupopt.problems.esipp.blocks.prices import NODE_PRICE_OTHER, NODE_PRICE_DELTA, NODE_PRICE_LAMBDA
 
 # *****************************************************************************
 # *****************************************************************************
@@ -29,8 +30,16 @@ class TestESIPPProblem:
     # *************************************************************************
     # *************************************************************************
 
-    @pytest.mark.parametrize("use_prices_block", [True, False])
-    def test_problem_increasing_imp_prices(self, use_prices_block):
+    @pytest.mark.parametrize(
+        "use_prices_block, node_price_model", 
+        [(True, NODE_PRICE_OTHER),
+         (True, NODE_PRICE_DELTA),
+         (True, NODE_PRICE_LAMBDA),
+         (False, NODE_PRICE_OTHER),
+         (False, NODE_PRICE_DELTA),
+         (False, NODE_PRICE_LAMBDA)]
+        )
+    def test_problem_increasing_imp_prices(self, use_prices_block, node_price_model):
         
         # assessment
         q = 0
@@ -48,10 +57,12 @@ class TestESIPPProblem:
 
         # import node
         node_IMP = 'I'
+        prices = [1.0, 2.0]
+        volumes = [0.5, None] if node_price_model == NODE_PRICE_OTHER else [0.5, 1e5]
         mynet.add_import_node(
             node_key=node_IMP,
             prices={
-                qpk: ResourcePrice(prices=[1.0, 2.0], volumes=[0.5, None])
+                qpk: ResourcePrice(prices=prices, volumes=volumes)
                 for qpk in tf.qpk()
             },
         )
@@ -86,10 +97,19 @@ class TestESIPPProblem:
             mandatory_arcs=[],
             max_number_parallel_arcs={},
             simplify_problem=False,
-            use_prices_block=use_prices_block
+            use_prices_block=use_prices_block,
+            node_price_model=node_price_model
         )
 
         assert not ipp.has_peak_total_assessments()
+        
+        print('hey')
+        print(ipp.results["Problem"][0])
+        # if node_price_model != NODE_PRICE_OTHER:
+        #     ipp.instance.pprint()
+        # assert ipp.results["Problem"][0]["Number of constraints"] == 24
+        # assert ipp.results["Problem"][0]["Number of variables"] == 22
+        # assert ipp.results["Problem"][0]["Number of nonzeros"] == 49
 
         # *********************************************************************
         # *********************************************************************
@@ -130,8 +150,16 @@ class TestESIPPProblem:
     # *************************************************************************
     # *************************************************************************
 
-    @pytest.mark.parametrize("use_prices_block", [True, False])
-    def test_problem_decreasing_imp_prices(self, use_prices_block):
+    @pytest.mark.parametrize(
+        "use_prices_block, node_price_model", 
+        [(True, NODE_PRICE_OTHER),
+         (True, NODE_PRICE_DELTA),
+         (True, NODE_PRICE_LAMBDA),
+         (False, NODE_PRICE_OTHER),
+         (False, NODE_PRICE_DELTA),
+         (False, NODE_PRICE_LAMBDA)]
+        )
+    def test_problem_decreasing_imp_prices(self, use_prices_block, node_price_model):
         
         # assessment
         q = 0
@@ -187,7 +215,8 @@ class TestESIPPProblem:
             mandatory_arcs=[],
             max_number_parallel_arcs={},
             simplify_problem=False,
-            use_prices_block=use_prices_block
+            use_prices_block=use_prices_block,
+            node_price_model=node_price_model
         )
 
         assert not ipp.has_peak_total_assessments()
@@ -231,8 +260,16 @@ class TestESIPPProblem:
     # *************************************************************************
     # *************************************************************************
 
-    @pytest.mark.parametrize("use_prices_block", [True, False])
-    def test_problem_decreasing_imp_prices_infinite_capacity(self, use_prices_block):
+    @pytest.mark.parametrize(
+        "use_prices_block, node_price_model", 
+        [(True, NODE_PRICE_OTHER),
+         (True, NODE_PRICE_DELTA),
+         (True, NODE_PRICE_LAMBDA),
+         (False, NODE_PRICE_OTHER),
+         (False, NODE_PRICE_DELTA),
+         (False, NODE_PRICE_LAMBDA)]
+        )
+    def test_problem_decreasing_imp_prices_infinite_capacity(self, use_prices_block, node_price_model):
         
         # assessment
         q = 0
@@ -299,8 +336,16 @@ class TestESIPPProblem:
     # *************************************************************************
     # *************************************************************************
 
-    @pytest.mark.parametrize("use_prices_block", [True, False])
-    def test_problem_decreasing_exp_prices(self, use_prices_block):
+    @pytest.mark.parametrize(
+        "use_prices_block, node_price_model", 
+        [(True, NODE_PRICE_OTHER),
+         (True, NODE_PRICE_DELTA),
+         (True, NODE_PRICE_LAMBDA),
+         (False, NODE_PRICE_OTHER),
+         (False, NODE_PRICE_DELTA),
+         (False, NODE_PRICE_LAMBDA)]
+        )
+    def test_problem_decreasing_exp_prices(self, use_prices_block, node_price_model):
         # assessment
         q = 0
         # time
@@ -321,10 +366,12 @@ class TestESIPPProblem:
 
         # import node
         node_EXP = generate_pseudo_unique_key(mynet.nodes())
+        prices = [2.0, 1.0]
+        volumes = [0.5, None] if node_price_model == NODE_PRICE_OTHER else [0.5, 1e5]
         mynet.add_export_node(
             node_key=node_EXP,
             prices={
-                (q, p, k): ResourcePrice(prices=[2.0, 1.0], volumes=[0.5, None])
+                (q, p, k): ResourcePrice(prices=prices, volumes=volumes)
                 for p in range(number_periods)
                 for k in range(number_intervals)
             },
@@ -360,6 +407,8 @@ class TestESIPPProblem:
             mandatory_arcs=[],
             max_number_parallel_arcs={},
             simplify_problem=False,
+            use_prices_block=use_prices_block,
+            node_price_model=node_price_model
         )
 
         assert not ipp.has_peak_total_assessments()
@@ -403,8 +452,16 @@ class TestESIPPProblem:
     # *************************************************************************
     # *************************************************************************
 
-    @pytest.mark.parametrize("use_prices_block", [True, False])
-    def test_problem_increasing_exp_prices(self, use_prices_block):
+    @pytest.mark.parametrize(
+        "use_prices_block, node_price_model", 
+        [(True, NODE_PRICE_OTHER),
+         (True, NODE_PRICE_DELTA),
+         (True, NODE_PRICE_LAMBDA),
+         (False, NODE_PRICE_OTHER),
+         (False, NODE_PRICE_DELTA),
+         (False, NODE_PRICE_LAMBDA)]
+        )
+    def test_problem_increasing_exp_prices(self, use_prices_block, node_price_model):
         # assessment
         q = 0
         # time
@@ -464,6 +521,8 @@ class TestESIPPProblem:
             mandatory_arcs=[],
             max_number_parallel_arcs={},
             simplify_problem=False,
+            use_prices_block=use_prices_block,
+            node_price_model=node_price_model
         )
         
         assert not ipp.has_peak_total_assessments()
@@ -507,8 +566,16 @@ class TestESIPPProblem:
     # *************************************************************************
     # *************************************************************************
 
-    @pytest.mark.parametrize("use_prices_block", [True, False])
-    def test_problem_increasing_exp_prices_infinite_capacity(self, use_prices_block):
+    @pytest.mark.parametrize(
+        "use_prices_block, node_price_model", 
+        [(True, NODE_PRICE_OTHER),
+         (True, NODE_PRICE_DELTA),
+         (True, NODE_PRICE_LAMBDA),
+         (False, NODE_PRICE_OTHER),
+         (False, NODE_PRICE_DELTA),
+         (False, NODE_PRICE_LAMBDA)]
+        )
+    def test_problem_increasing_exp_prices_infinite_capacity(self, use_prices_block, node_price_model):
         # assessment
         q = 0
         # time
@@ -571,7 +638,8 @@ class TestESIPPProblem:
                 mandatory_arcs=[],
                 max_number_parallel_arcs={},
                 simplify_problem=False,
-                use_prices_block=use_prices_block
+                use_prices_block=use_prices_block,
+                node_price_model=node_price_model
             )
         except Exception:
             error_raised = True
@@ -580,8 +648,16 @@ class TestESIPPProblem:
     # *************************************************************************
     # *************************************************************************
 
-    @pytest.mark.parametrize("use_prices_block", [True, False])
-    def test_problem_increasing_imp_decreasing_exp_prices(self, use_prices_block):
+    @pytest.mark.parametrize(
+        "use_prices_block, node_price_model", 
+        [(True, NODE_PRICE_OTHER),
+         (True, NODE_PRICE_DELTA),
+         (True, NODE_PRICE_LAMBDA),
+         (False, NODE_PRICE_OTHER),
+         (False, NODE_PRICE_DELTA),
+         (False, NODE_PRICE_LAMBDA)]
+        )
+    def test_problem_increasing_imp_decreasing_exp_prices(self, use_prices_block, node_price_model):
         # scenario
         q = 0
         # time
@@ -602,10 +678,12 @@ class TestESIPPProblem:
 
         # import node
         node_IMP = 'I'
+        prices = [1.0, 2.0]
+        volumes = [0.5, None] if node_price_model == NODE_PRICE_OTHER else [0.5, 1e5]
         mynet.add_import_node(
             node_key=node_IMP,
             prices={
-                (q, p, k): ResourcePrice(prices=[1.0, 2.0], volumes=[0.5, None])
+                (q, p, k): ResourcePrice(prices=prices, volumes=volumes)
                 for p in range(number_periods)
                 for k in range(number_intervals)
             },
@@ -613,10 +691,12 @@ class TestESIPPProblem:
 
         # export node
         node_EXP = generate_pseudo_unique_key(mynet.nodes())
+        prices = [2.0, 1.0]
+        volumes = [0.5, None] if node_price_model == NODE_PRICE_OTHER else [0.5, 1e5]
         mynet.add_export_node(
             node_key=node_EXP,
             prices={
-                (q, p, k): ResourcePrice(prices=[2.0, 1.0], volumes=[0.5, None])
+                (q, p, k): ResourcePrice(prices=prices, volumes=volumes)
                 for p in range(number_periods)
                 for k in range(number_intervals)
             },
@@ -669,7 +749,8 @@ class TestESIPPProblem:
             max_number_parallel_arcs={},
             simplify_problem=False,
             # discount_rates={0: (0.0,)},
-            use_prices_block=use_prices_block
+            use_prices_block=use_prices_block,
+            node_price_model=node_price_model
         )
 
         assert not ipp.has_peak_total_assessments()
@@ -743,8 +824,16 @@ class TestESIPPProblem:
     # *************************************************************************
     # *************************************************************************
     
-    @pytest.mark.parametrize("use_prices_block", [True, False])
-    def test_direct_imp_exp_network_higher_exp_prices(self, use_prices_block):
+    @pytest.mark.parametrize(
+        "use_prices_block, node_price_model", 
+        [(True, NODE_PRICE_OTHER),
+         (True, NODE_PRICE_DELTA),
+         (True, NODE_PRICE_LAMBDA),
+         (False, NODE_PRICE_OTHER),
+         (False, NODE_PRICE_DELTA),
+         (False, NODE_PRICE_LAMBDA)]
+        )
+    def test_direct_imp_exp_network_higher_exp_prices(self, use_prices_block, node_price_model):
         
         # time frame
         q = 0
@@ -764,7 +853,7 @@ class TestESIPPProblem:
         imp_prices = {
             qpk: ResourcePrice(
                 prices=0.5,
-                volumes=None,
+                volumes=None if node_price_model == NODE_PRICE_OTHER else 1e4,
             )
             for qpk in tf.qpk()
             }
@@ -778,7 +867,7 @@ class TestESIPPProblem:
         exp_prices = {
             qpk: ResourcePrice(
                 prices=1.5,
-                volumes=None,
+                volumes=None if node_price_model == NODE_PRICE_OTHER else 1e4,
             )
             for qpk in tf.qpk()
             }
@@ -815,7 +904,8 @@ class TestESIPPProblem:
             static_losses_mode=InfrastructurePlanningProblem.STATIC_LOSS_MODE_DEP,
             mandatory_arcs=[],
             max_number_parallel_arcs={},
-            use_prices_block=use_prices_block
+            use_prices_block=use_prices_block,
+            node_price_model=node_price_model
         )
     
         # export prices are higher: it makes sense to install the arc since the
diff --git a/tests/test_esipp_problem.py b/tests/test_esipp_problem.py
index a1072d36aba25e99bf9a862065ed77a30ef57fae..d1005617716abd3a8fd761b88d469d462dfdd291 100644
--- a/tests/test_esipp_problem.py
+++ b/tests/test_esipp_problem.py
@@ -824,8 +824,17 @@ class TestESIPPProblem:
     # problem with symmetrical nodes and one undirected arc, irregular steps
     # same problem as the previous one, except with interface variables
     # problem with two symmetrical nodes and one undirected arc, w/ simple sos1
-    # TODO: test SOS 
-    def test_isolated_undirected_network(self):
+    
+    @pytest.mark.parametrize(
+        "solver, use_sos_arcs, arc_sos_weight_key", 
+        [('scip', True, InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_NONE),
+         # ('scip', True, InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_CAP),
+         # ('scip', True, InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_COST),
+         # ('scip', True, InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_SPEC_CAP),
+         # ('scip', True, InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_SPEC_COST),
+         ('scip', False, None)]
+        )
+    def test_isolated_undirected_network(self, solver, use_sos_arcs, arc_sos_weight_key):
         
         q = 0
         tf = EconomicTimeFrame(
@@ -873,6 +882,8 @@ class TestESIPPProblem:
         solver = 'scip'
         ipp = build_solve_ipp(
             solver=solver,
+            use_sos_arcs=use_sos_arcs,
+            arc_sos_weight_key=arc_sos_weight_key,
             solver_options={},
             perform_analysis=False,
             plot_results=False,  # True,
@@ -884,7 +895,11 @@ class TestESIPPProblem:
             max_number_parallel_arcs={}
         )
         
-        assert len(ipp.instance.constr_arc_sos1) == 0
+        if use_sos_arcs:
+            assert len(ipp.instance.constr_arc_sos1) != 0
+        else:
+            assert len(ipp.instance.constr_arc_sos1) == 0
+            
         assert ipp.has_peak_total_assessments() # TODO: make sure this is true
         # assert ipp.results["Problem"][0]["Number of constraints"] == 34
         # assert ipp.results["Problem"][0]["Number of variables"] == 28
@@ -2512,7 +2527,6 @@ class TestESIPPProblem:
     # *************************************************************************
     # *************************************************************************
     
-    # TODO: test arc groups with sos
     @pytest.mark.parametrize(
         "solver, use_sos_arc_groups, arc_sos_weight_key", 
         [('scip', True, InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_CAP),
@@ -7415,7 +7429,6 @@ class TestESIPPProblem:
     # TODO: test non-simplifiable problems with time varying prices on select assessments
     # TODO: test non-simplifiable problems with volume varying prices on select assessments
     
-        
     # *************************************************************************
     # *************************************************************************
     
@@ -7525,9 +7538,6 @@ class TestESIPPProblem:
             simplify_problem=True,
         )
         assert ipp.has_peak_total_assessments()
-        # assert ipp.results["Problem"][0]["Number of constraints"] == 61
-        # assert ipp.results["Problem"][0]["Number of variables"] == 53 
-        # assert ipp.results["Problem"][0]["Number of nonzeros"] == 143
         
         # *********************************************************************
         # *********************************************************************
@@ -7669,9 +7679,6 @@ class TestESIPPProblem:
             simplify_problem=True,
         )
         assert ipp.has_peak_total_assessments()
-        # assert ipp.results["Problem"][0]["Number of constraints"] == 61
-        # assert ipp.results["Problem"][0]["Number of variables"] == 53 
-        # assert ipp.results["Problem"][0]["Number of nonzeros"] == 143 # 
         
         # *********************************************************************
         # *********************************************************************