diff --git a/src/topupopt/problems/esipp/blocks/prices.py b/src/topupopt/problems/esipp/blocks/prices.py
index 48aae65e479ff7737c9096bc039c55ad87bc0582..8fc17ca0d1f5df77e235b10251cdd385d81422b3 100644
--- a/src/topupopt/problems/esipp/blocks/prices.py
+++ b/src/topupopt/problems/esipp/blocks/prices.py
@@ -205,7 +205,7 @@ def price_block_other(
     # if not convex_price_function:
         
     #     # delta variables
-    #     model.var_delta_glqpks = pyo.Var(
+    #     model.var_active_segment_glqpks = pyo.Var(
     #         model.set_GLQPKS, within=pyo.Binary
     #     )
         
@@ -214,7 +214,7 @@ def price_block_other(
     #         return (
     #             m.var_if_glqpks[(g,l,q,p,k,s)] <= 
     #             m.param_v_max_glqpks[(g,l,q,p,k,s)]*
-    #             m.var_delta_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_imp = pyo.Constraint(
     #         model.set_GLQPKS_imp, rule=rule_constr_empty_segment_if_delta_zero_imp
@@ -225,7 +225,7 @@ def price_block_other(
     #         return (
     #             m.var_ef_glqpks[(g,l,q,p,k,s)] <= 
     #             m.param_v_max_glqpks[(g,l,q,p,k,s)]*
-    #             m.var_delta_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_exp = pyo.Constraint(
     #         model.set_GLQPKS_exp, rule=rule_constr_empty_segment_if_delta_zero_exp
@@ -236,8 +236,8 @@ def price_block_other(
     #         if s == len(m.set_S)-1:
     #             return pyo.Constraint.Skip
     #         return (
-    #             m.var_delta_glqpks[(g,l,q,p,k,s)] >= 
-    #             m.var_delta_glqpks[(g,l,q,p,k,s+1)]
+    #             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
@@ -247,8 +247,8 @@ def price_block_other(
     #         if s == len(m.set_S)-1:
     #             return pyo.Constraint.Skip
     #         return (
-    #             1-m.var_delta_glqpks[(g,l,q,p,k,s)] >= 
-    #             m.var_delta_glqpks[(g,l,q,p,k,s+1)]
+    #             1-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_next_zeros = pyo.Constraint(
     #         model.set_GLQPKS, rule=rule_constr_delta_next_zeros
@@ -256,20 +256,261 @@ def price_block_other(
 
     # *************************************************************************
     # *************************************************************************
+    
+    
+# *****************************************************************************
+# *****************************************************************************
+
+# def price_other2(
+#     model: pyo.AbstractModel,
+#     convex_price_function: bool = False,
+#     enable_default_values: bool = True,
+#     enable_validation: bool = True,
+#     enable_initialisation: bool = True
+#     ):
+
+#     # set of price segments
+#     model.set_S = pyo.Set(model.set_GL_exp_imp, model.set_QPK)
+
+#     # 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)
+#     )
+
+#     def init_set_GLQPKS_exp(m):
+#         return (
+#             glqpks for glqpks in m.set_GLQPKS if glqpks[1] in m.set_L_exp[glqpks[0]]
+#         )
+
+#     model.set_GLQPKS_exp = pyo.Set(
+#         dimen=6, initialize=(init_set_GLQPKS_exp if enable_initialisation else None)
+#     )
+
+#     def init_set_GLQPKS_imp(m):
+#         return (
+#             glqpks for glqpks in m.set_GLQPKS if glqpks[1] in m.set_L_imp[glqpks[0]]
+#         )
+
+#     model.set_GLQPKS_imp = pyo.Set(
+#         dimen=6, initialize=(init_set_GLQPKS_imp if enable_initialisation else None)
+#     )
+
+#     # *************************************************************************
+#     # *************************************************************************
+
+#     # parameters
+
+#     # resource prices
+
+#     model.param_p_glqpks = pyo.Param(model.set_GLQPKS, within=pyo.NonNegativeReals)
+
+#     # maximum resource volumes for each prices
+
+#     model.param_v_max_glqpks = pyo.Param(
+#         model.set_GLQPKS, 
+#         within=pyo.NonNegativeReals
+#         )
+
+#     # *************************************************************************
+#     # *************************************************************************
+
+#     # variables
+
+#     # *************************************************************************
+#     # *************************************************************************
+
+#     # exported flow
+
+#     # TODO: validate the bounds by ensuring inf. cap. only exists in last segm.
+
+#     def bounds_var_ef_glqpks(m, 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_ef_glqpks = pyo.Var(
+#         model.set_GLQPKS_exp, within=pyo.NonNegativeReals, bounds=bounds_var_ef_glqpks
+#     )
+
+#     # imported flow
+
+#     def bounds_var_if_glqpks(m, 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_if_glqpks = pyo.Var(
+#         model.set_GLQPKS_imp, within=pyo.NonNegativeReals, bounds=bounds_var_if_glqpks
+#     )
+
+#     # *************************************************************************
+#     # *************************************************************************
+
+#     # exported flow revenue
+#     def rule_constr_exp_flow_revenue(m, g, l, q, p, k):
+#         return (
+#             sum(
+#                 m.var_ef_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_exp_flow_revenue = pyo.Constraint(
+#         model.set_GL_exp, model.set_QPK, rule=rule_constr_exp_flow_revenue
+#     )
+
+#     # imported flow cost
+#     def rule_constr_imp_flow_cost(m, g, l, q, p, k):
+#         return (
+#             sum(
+#                 m.var_if_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)]
+#         )
+
+#     model.constr_imp_flow_cost = pyo.Constraint(
+#         model.set_GL_imp, model.set_QPK, rule=rule_constr_imp_flow_cost
+#     )
+
+#     # exported flows
+#     def rule_constr_exp_flows(m, g, l, q, p, k):
+#         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_ef_glqpks[(g, l, q, p, k, s)] for s in m.set_S[(g, l, q, p, k)])
+
+#     model.constr_exp_flows = pyo.Constraint(
+#         model.set_GL_exp, model.set_QPK, rule=rule_constr_exp_flows
+#     )
+
+#     # imported flows
+#     def rule_constr_imp_flows(m, g, l, q, p, k):
+#         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_if_glqpks[(g, l, q, p, k, s)] for s in m.set_S[(g, l, q, p, k)])
+
+#     model.constr_imp_flows = pyo.Constraint(
+#         model.set_GL_imp, model.set_QPK, rule=rule_constr_imp_flows
+#     )
+
+#     # *************************************************************************
+#     # *************************************************************************
+    
+#     # non-convex price functions
+    
+#     if not convex_price_function:
+        
+#         # 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_imp(m, 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)]
+#                 )
+#         model.constr_empty_segment_if_delta_zero_imp = pyo.Constraint(
+#             model.set_GLQPKS_imp, rule=rule_constr_empty_segment_if_delta_zero_imp
+#             )
+            
+#         # segments must be empty if the respective delta variable is zero
+#         def rule_constr_empty_segment_if_delta_zero_exp(m, g, l, q, p, k, s):
+#             return (
+#                 m.var_ef_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_exp = pyo.Constraint(
+#             model.set_GLQPKS_exp, rule=rule_constr_empty_segment_if_delta_zero_exp
+#             )
+        
+#         # 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:
+#                 # last segment, 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:
+#                 # last segment, skip
+#                 return pyo.Constraint.Skip
+#             if (g,l) in m.set_GL_imp:
+#                 return (
+#                     m.var_if_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)]
+#                     )
+#             else:
+#                 return (
+#                     m.var_ef_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
+#             )
 
 # *****************************************************************************
 # *****************************************************************************
 
 def price_other(
     model: pyo.AbstractModel,
-    convex_price_function: bool = False,
+    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_GL_exp_imp, model.set_QPK)
+    model.set_S = pyo.Set(model.set_GLQPK)
 
     # set of GLQKS tuples
     def init_set_GLQPKS(m):
@@ -285,24 +526,6 @@ def price_other(
         dimen=6, initialize=(init_set_GLQPKS if enable_initialisation else None)
     )
 
-    def init_set_GLQPKS_exp(m):
-        return (
-            glqpks for glqpks in m.set_GLQPKS if glqpks[1] in m.set_L_exp[glqpks[0]]
-        )
-
-    model.set_GLQPKS_exp = pyo.Set(
-        dimen=6, initialize=(init_set_GLQPKS_exp if enable_initialisation else None)
-    )
-
-    def init_set_GLQPKS_imp(m):
-        return (
-            glqpks for glqpks in m.set_GLQPKS if glqpks[1] in m.set_L_imp[glqpks[0]]
-        )
-
-    model.set_GLQPKS_imp = pyo.Set(
-        dimen=6, initialize=(init_set_GLQPKS_imp if enable_initialisation else None)
-    )
-
     # *************************************************************************
     # *************************************************************************
 
@@ -311,6 +534,13 @@ def price_other(
     # 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
 
@@ -327,172 +557,129 @@ def price_other(
     # *************************************************************************
     # *************************************************************************
 
-    # exported flow
-
-    # TODO: validate the bounds by ensuring inf. cap. only exists in last segm.
-
-    def bounds_var_ef_glqpks(m, g, l, q, p, k, s):
+    # import and export flows
+    def bounds_var_trans_flows_glqpks(m, 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_ef_glqpks = pyo.Var(
-        model.set_GLQPKS_exp, within=pyo.NonNegativeReals, bounds=bounds_var_ef_glqpks
-    )
-
-    # imported flow
-
-    def bounds_var_if_glqpks(m, 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_if_glqpks = pyo.Var(
-        model.set_GLQPKS_imp, within=pyo.NonNegativeReals, bounds=bounds_var_if_glqpks
+    model.var_trans_flows_glqpks = pyo.Var(
+        model.set_GLQPKS, within=pyo.NonNegativeReals, bounds=bounds_var_trans_flows_glqpks
     )
 
     # *************************************************************************
     # *************************************************************************
 
-    # exported flow revenue
-    def rule_constr_exp_flow_revenue(m, g, l, q, p, k):
-        return (
-            sum(
-                m.var_ef_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)]
+    # 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)]
             )
-            == m.var_efr_glqpk[(g, l, q, p, k)]
-        )
-
-    model.constr_exp_flow_revenue = pyo.Constraint(
-        model.set_GL_exp, model.set_QPK, rule=rule_constr_exp_flow_revenue
-    )
-
-    # imported flow cost
-    def rule_constr_imp_flow_cost(m, g, l, q, p, k):
-        return (
-            sum(
-                m.var_if_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)]
+        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)]
             )
-            == m.var_ifc_glqpk[(g, l, q, p, k)]
-        )
-
-    model.constr_imp_flow_cost = pyo.Constraint(
-        model.set_GL_imp, model.set_QPK, rule=rule_constr_imp_flow_cost
+    model.constr_trans_monetary_flows = pyo.Constraint(
+        model.set_GLQPK, rule=rule_constr_trans_monetary_flows
     )
 
-    # exported flows
-    def rule_constr_exp_flows(m, g, l, q, p, k):
-        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_ef_glqpks[(g, l, q, p, k, s)] for s in m.set_S[(g, l, q, p, k)])
-
-    model.constr_exp_flows = pyo.Constraint(
-        model.set_GL_exp, model.set_QPK, rule=rule_constr_exp_flows
-    )
-
-    # imported flows
-    def rule_constr_imp_flows(m, g, l, q, p, k):
-        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_if_glqpks[(g, l, q, p, k, s)] for s in m.set_S[(g, l, q, p, k)])
+    # 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_imp_flows = pyo.Constraint(
-        model.set_GL_imp, model.set_QPK, rule=rule_constr_imp_flows
+    model.constr_trans_flows = pyo.Constraint(
+        model.set_GLQPK, rule=rule_constr_trans_flows
     )
 
     # *************************************************************************
     # *************************************************************************
     
     # non-convex price functions
-    
-    if not convex_price_function:
-        
-        # delta variables
-        model.var_delta_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_imp(m, 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_delta_glqpks[(g,l,q,p,k,s)]
-                )
-        model.constr_empty_segment_if_delta_zero_imp = pyo.Constraint(
-            model.set_GLQPKS_imp, rule=rule_constr_empty_segment_if_delta_zero_imp
-            )
-            
-        # segments must be empty if the respective delta variable is zero
-        def rule_constr_empty_segment_if_delta_zero_exp(m, g, l, q, p, k, s):
-            return (
-                m.var_ef_glqpks[(g,l,q,p,k,s)] <= 
-                m.param_v_max_glqpks[(g,l,q,p,k,s)]*
-                m.var_delta_glqpks[(g,l,q,p,k,s)]
-                )
-        model.constr_empty_segment_if_delta_zero_exp = pyo.Constraint(
-            model.set_GLQPKS_exp, rule=rule_constr_empty_segment_if_delta_zero_exp
+    # 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)]
             )
-        
-        # 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:
-                # last segment, skip
-                return pyo.Constraint.Skip
-            return (
-                m.var_delta_glqpks[(g,l,q,p,k,s)] >= 
-                m.var_delta_glqpks[(g,l,q,p,k,s+1)]
-                )
-        model.constr_delta_summing_logic = pyo.Constraint(
-            model.set_GLQPKS, rule=rule_constr_delta_summing_logic
+    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)]
             )
-        
-        # 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:
-                # last segment, skip
-                return pyo.Constraint.Skip
-            if (g,l) in m.set_GL_imp:
-                return (
-                    m.var_if_glqpks[(g,l,q,p,k,s)] >= 
-                    m.var_delta_glqpks[(g,l,q,p,k,s+1)]*
-                    m.param_v_max_glqpks[(g,l,q,p,k,s)]
-                    )
-            else:
-                return (
-                    m.var_ef_glqpks[(g,l,q,p,k,s)] >= 
-                    m.var_delta_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_delta_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_delta_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
+    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
+        )
 
 # *****************************************************************************
 # *****************************************************************************
diff --git a/src/topupopt/problems/esipp/problem.py b/src/topupopt/problems/esipp/problem.py
index 5b2a5707492400360d4b7c66ff1df11c319b87f9..3659ad67caea582d20dab8d2251aa3d023987d71 100644
--- a/src/topupopt/problems/esipp/problem.py
+++ b/src/topupopt/problems/esipp/problem.py
@@ -2549,6 +2549,17 @@ class InfrastructurePlanningProblem(EnergySystem):
             for s in set_S[(g, l, q, p, k)]
         }
 
+        # price function convexity
+
+        param_price_function_is_convex = {
+            (g, l, q, p, k): (
+                self.networks[g].nodes[l][Network.KEY_NODE_PRICES][(q, p, k)].price_monotonically_increasing_with_volume()
+                if l in set_L_imp[g] else
+                self.networks[g].nodes[l][Network.KEY_NODE_PRICES][(q, p, k)].price_monotonically_decreasing_with_volume()
+                )
+            for (g, l, q, p, k) in set_S
+        }
+
         # maximum resource volume per segment (infinity is the default)
 
         param_v_max_glqpks = {
@@ -3451,6 +3462,7 @@ class InfrastructurePlanningProblem(EnergySystem):
                 "param_c_df_qp": param_c_df_qp,
                 "param_c_time_qpk": param_c_time_qpk,
                 "param_p_glqpks": param_p_glqpks,
+                "param_price_function_is_convex": param_price_function_is_convex,
                 "param_v_max_glqpks": param_v_max_glqpks,
                 # *****************************************************************
                 # converters
diff --git a/src/topupopt/problems/esipp/utils.py b/src/topupopt/problems/esipp/utils.py
index 94c3f790cc9049638b683fd2dc4443cabdc0d323..4047696c1311ed28b21307510000cd5c8638cbf9 100644
--- a/src/topupopt/problems/esipp/utils.py
+++ b/src/topupopt/problems/esipp/utils.py
@@ -99,7 +99,7 @@ def statistics(ipp: InfrastructurePlanningProblem,
     imports_qpk = {
         qpk: pyo.value( 
             sum(
-                ipp.instance.var_if_glqpks[(g,l_imp,*qpk, s)]
+                ipp.instance.var_trans_flows_glqpks[(g,l_imp,*qpk, s)]
                 for g, l_imp in import_node_keys
                 # for g in ipp.networks
                 # for l_imp in ipp.networks[g].import_nodes
@@ -114,7 +114,7 @@ def statistics(ipp: InfrastructurePlanningProblem,
     exports_qpk = {
         qpk: pyo.value(
             sum(
-                ipp.instance.var_ef_glqpks[(g,l_exp,*qpk, s)]
+                ipp.instance.var_trans_flows_glqpks[(g,l_exp,*qpk, s)]
                 for g, l_exp in export_node_keys
                 # for g in ipp.networks
                 # for l_exp in ipp.networks[g].export_nodes
diff --git a/tests/test_esipp_prices.py b/tests/test_esipp_prices.py
index 31e4c55d7b5638d1a10fe0cf5896545cf1656f72..783304755c4d764ac00cd70e0ac0da3deb620f71 100644
--- a/tests/test_esipp_prices.py
+++ b/tests/test_esipp_prices.py
@@ -629,8 +629,6 @@ class TestESIPPProblem:
         # the objective function should be -7.5
         assert math.isclose(pyo.value(ipp.instance.obj_f), -2.0, abs_tol=1e-3)
         
-    # TODO: trigger errors with infinite capacity segments in non-convex functions
-        
     # *************************************************************************
     # *************************************************************************