Added support for convex and nonconvex price functions using the lambda...

Added support for convex and nonconvex price functions using the lambda formulation, now also without blocks. Tests were slightly adjusted too.

src/topupopt/problems/esipp/blocks/lambda_.py

@@ -49,29 +49,29 @@ def price_lambda_block(
         # *********************************************************************
         
         # set of points
-        b.set_P = pyo.Set(initialize=[i for i in range(len(b.set_S)+1)])
+        b.set_O = pyo.Set(initialize=[i for i in range(len(b.set_S)+1)])
         
         # set of x points
-        def init_param_x_p(b, p):
-            return 0 if p == 0 else sum(b.param_v_max_s[i] for i in range(p))
-        b.param_x_p = pyo.Param(
-            b.set_P, 
+        def init_param_x_o(b, o):
+            return 0 if o == 0 else sum(b.param_v_max_s[i] for i in range(o))
+        b.param_x_o = pyo.Param(
+            b.set_O, 
             within=pyo.NonNegativeReals,
-            initialize=init_param_x_p
+            initialize=init_param_x_o
             )
         
         # set of y points
-        def init_param_y_p(b, p):
-            return 0 if p == 0 else (b.param_x_p[p]-b.param_x_p[p-1])*b.param_p_s[p-1]+b.param_y_p[p-1]
+        def init_param_y_o(b, o):
+            return 0 if o == 0 else (b.param_x_o[o]-b.param_x_o[o-1])*b.param_p_s[o-1]+b.param_y_o[o-1]
         # sum(b.param_p_s[i] for i in range(p+1))
-        b.param_y_p = pyo.Param(
-            b.set_P, 
+        b.param_y_o = pyo.Param(
+            b.set_O, 
             within=pyo.NonNegativeReals,
-            initialize=init_param_y_p
+            initialize=init_param_y_o
             )
     
         # interpoint weights
-        b.var_weights_p = pyo.Var(b.set_P, within=pyo.UnitInterval)
+        b.var_weights_o = pyo.Var(b.set_O, within=pyo.UnitInterval)
         
         # TODO: consider replacing None in the bounds with inf in the parameter
         
@@ -87,7 +87,7 @@ def price_lambda_block(
         # *********************************************************************
         
         def rule_constr_stick_to_line(b):
-            return sum(b.var_weights_p[p] for p in b.set_P) == 1
+            return sum(b.var_weights_o[p] for p in b.set_O) == 1
         b.constr_stick_to_line = pyo.Constraint(rule=rule_constr_stick_to_line)
         
         # *********************************************************************
@@ -97,12 +97,12 @@ def price_lambda_block(
         def rule_constr_y_equation(b):
             if (g,l) in b.parent_block().set_GL_imp:
                 return (
-                    sum(b.var_weights_p[p]*b.param_y_p[p] for p in b.set_P)
+                    sum(b.var_weights_o[p]*b.param_y_o[p] for p in b.set_O)
                     == b.parent_block().var_ifc_glqpk[(g,l,q,p,k)]
                 )
             else:
                 return (
-                    sum(b.var_weights_p[p]*b.param_y_p[p] for p in b.set_P)
+                    sum(b.var_weights_o[p]*b.param_y_o[p] for p in b.set_O)
                     == b.parent_block().var_efr_glqpk[(g,l,q,p,k)]
                 )
         b.constr_y_equation = pyo.Constraint(rule=rule_constr_y_equation)
@@ -111,12 +111,12 @@ def price_lambda_block(
         def rule_constr_x_equation(b):
             if (g,l) in b.parent_block().set_GL_imp:
                 return (
-                    sum(b.var_weights_p[p]*b.param_x_p[p] for p in b.set_P)
+                    sum(b.var_weights_o[p]*b.param_x_o[p] for p in b.set_O)
                     == b.var_trans_flows
                 )
             else:
                 return (
-                    sum(b.var_weights_p[p]*b.param_x_p[p] for p in b.set_P)
+                    sum(b.var_weights_o[p]*b.param_x_o[p] for p in b.set_O)
                     == b.var_trans_flows
                 )
         b.constr_x_equation = pyo.Constraint(rule=rule_constr_x_equation)
@@ -148,7 +148,7 @@ def price_lambda_block(
             
             # declare SOS2
             def rule_constr_sos2_weights(b):
-                return [b.var_weights_p[p] for p in b.set_P]            
+                return [b.var_weights_o[p] for p in b.set_O]            
             b.constr_sos2_weights = pyo.SOSConstraint(rule=rule_constr_sos2_weights, sos=2)
         
         # *********************************************************************
@@ -180,8 +180,6 @@ def price_lambda_no_block(
     enable_initialisation: bool = True
     ):
     
-    raise NotImplementedError
-    
     # auxiliary set for pyomo
     model.set_GLQPK = model.set_GL_exp_imp*model.set_QPK
 
@@ -208,18 +206,15 @@ def price_lambda_no_block(
     # 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,
@@ -234,55 +229,135 @@ def price_lambda_no_block(
     # *************************************************************************
     # *************************************************************************
     
-    # 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
+    # set of points
+    def init_set_O(m, g, l, q, p, k):
+        return [i for i in range(len(m.set_S[(g,l,q,p,k)])+1)]
+        # return 0 if o == 0 else sum(m.param_v_max_glqpks[i] for i in range(o))
+    model.set_O = pyo.Set(
+        model.set_GLQPK,
+        initialize=(init_set_O if enable_initialisation else None)
+        )
+    
+    
+    # set of GLQKS tuples
+    def init_set_GLQPKO(m):
+        return (
+            (g, l, q, p, k, o)
+            # 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_O
+            for o in m.set_O[(g, l, q, p, k)]
+        )
+    model.set_GLQPKO = pyo.Set(
+        dimen=6, initialize=(init_set_GLQPKO if enable_initialisation else None)
+    )
+    
+    # set of x points
+    def init_param_x_glqpko(m, g, l, q, p, k, o):
+        return 0 if o == 0 else sum(m.param_v_max_glqpks[(g,l,q,p,k,i)] for i in range(o))
+    model.param_x_glqpko = pyo.Param(
+        model.set_GLQPKO, 
+        within=pyo.NonNegativeReals,
+        initialize=init_param_x_glqpko
+        )
+    
+    # set of y points
+    def init_param_y_glqpko(m, g, l, q, p, k, o):
+        return (
+            0 
+            if o == 0 else 
+            (m.param_x_glqpko[(g,l,q,p,k,o)]
+             -m.param_x_glqpko[(g,l,q,p,k,o-1)])*
+            m.param_p_glqpks[(g,l,q,p,k,o-1)]+
+            m.param_y_glqpko[(g,l,q,p,k,o-1)]
+            )
+    model.param_y_glqpko = pyo.Param(
+        model.set_GLQPKO, 
+        within=pyo.NonNegativeReals,
+        initialize=init_param_y_glqpko
+        )
+
+    # interpoint weights
+    model.var_weights_glqpko = pyo.Var(
+        model.set_GLQPKO, 
+        within=pyo.UnitInterval
+        )
+    
+    # TODO: consider replacing None in the bounds with inf in the parameter
+    
+    # transshipment flows        
+    def bounds_var_trans_flows_glqpk(m, g, l, q, p, k):
+        try:
+            return (0, sum(m.param_v_max_glqpks[(g, l, q, p, k, s)]
+                           for s in m.set_S[(g, l, q, p, k)]))
+        except Exception:
             return (0, None)
-    model.var_trans_flows_glqpks = pyo.Var(
-        model.set_GLQPKS, within=pyo.NonNegativeReals, bounds=bounds_var_trans_flows_glqpks
+    model.var_trans_flows_glqpk = pyo.Var(
+        model.set_GLQPK, 
+        bounds=bounds_var_trans_flows_glqpk
+        )
+
+    # *************************************************************************
+    # *************************************************************************
+    
+    def rule_constr_stick_to_line(m, g, l, q, p, k):
+        return sum(m.var_weights_glqpko[(g,l,q,p,k,o)]
+                   for o in m.set_O[(g,l,q,p,k)]) == 1
+    model.constr_stick_to_line = pyo.Constraint(
+        model.set_GLQPK, 
+        rule=rule_constr_stick_to_line
         )
     
     # *************************************************************************
     # *************************************************************************
 
-    # import flow costs and export flow revenues
-    def rule_constr_trans_monetary_flows(m, g, l, q, p, k):
+    # import flow costs and export flow revenues (y equation)
+    def rule_constr_y_equation(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)]
-                )
+                sum(m.var_weights_glqpko[(g,l,q,p,k,o)]*
+                    m.param_y_glqpko[(g,l,q,p,k,o)] 
+                    for o in m.set_O[(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)]
-                )
+                sum(m.var_weights_glqpko[(g,l,q,p,k,o)]*
+                    m.param_y_glqpko[(g,l,q,p,k,o)] 
+                    for o in m.set_O[(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
+    model.constr_y_equation = pyo.Constraint(
+        model.set_GLQPK,
+        rule=rule_constr_y_equation
+        )
+
+    # imported and exported flows (x equation)
+    def rule_constr_x_equation(m, g, l, q, p, k):
+        if (g,l) in m.set_GL_imp:
+            return (
+                sum(m.var_weights_glqpko[(g,l,q,p,k,o)]*m.param_x_glqpko[(g,l,q,p,k,o)] for o in m.set_O[(g,l,q,p,k)])
+                == m.var_trans_flows_glqpk[(g,l,q,p,k)]
+            )
+        else:
+            return (
+                sum(m.var_weights_glqpko[(g,l,q,p,k,o)]*m.param_x_glqpko[(g,l,q,p,k,o)] for o in m.set_O[(g,l,q,p,k)])
+                == m.var_trans_flows_glqpk[(g,l,q,p,k)]
+            )
+    model.constr_x_equation = pyo.Constraint(
+        model.set_GLQPK,
+        rule=rule_constr_x_equation
         )
     
-    # imported and exported flows
-    def rule_constr_trans_flows(m, g, l, q, p, k):
+    # imported and exported flows (system equation)
+    def rule_constr_sys_equation(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)])
+            ) == m.var_trans_flows_glqpk[(g,l,q,p,k)]
         else:
             return sum(
                 m.var_v_glljqk[(g,l_star,l,j,q,k)]
@@ -290,74 +365,29 @@ def price_lambda_no_block(
                 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
+            ) == m.var_trans_flows_glqpk[(g,l,q,p,k)]
+    model.constr_sys_equation = pyo.Constraint(
+        model.set_GLQPK,
+        rule=rule_constr_sys_equation
         )
 
     # *************************************************************************
     # *************************************************************************
     
     # 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
-        )
+    # declare SOS2
+    def rule_constr_sos2_weights(m, g, l, q, p, k): 
+        if m.param_price_function_is_convex[(g,l,q,p,k)]:
+            return pyo.SOSConstraint.Skip
+        return [m.var_weights_glqpko[(g,l,q,p,k,o)] for o in m.set_O[(g,l,q,p,k)]]            
+    model.constr_sos2_weights = pyo.SOSConstraint(
+        model.set_GLQPK,
+        rule=rule_constr_sos2_weights, 
+        sos=2) 
 
-    # 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

tests/test_esipp_prices.py

@@ -163,8 +163,7 @@ class TestESIPPProblem:
         
     # *************************************************************************
     # *************************************************************************
-    # TODO: test a nonconvex price function in which only the first segment is used
-    # e.g. using a flow 0.5 in the following example
+    
     @pytest.mark.parametrize(
         "use_prices_block, node_price_model", 
         [(True, NODE_PRICE_OTHER),