diff --git a/src/topupopt/problems/esipp/blocks/lambda_.py b/src/topupopt/problems/esipp/blocks/lambda_.py
index 7cb66668206fc10acc26a501ad92a148f8ef9256..4a2d7223328d519f038aaddf4765eee2fc109cb6 100644
--- a/src/topupopt/problems/esipp/blocks/lambda_.py
+++ b/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,24 +206,21 @@ 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,
# default=math.inf
)
-
+
# *************************************************************************
# *************************************************************************
@@ -233,131 +228,166 @@ def price_lambda_no_block(
# *************************************************************************
# *************************************************************************
+
+ # 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
+ )
- # 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
+ # 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)]
- )
- == m.var_ifc_glqpk[(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)]
- )
- == m.var_efr_glqpk[(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
- def rule_constr_trans_flows(m, g, l, q, p, k):
+ # 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 (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)]
+ 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)])
+ for j in m.set_J[(g,l,l_star)] # only directed arcs
+ ) == m.var_trans_flows_glqpk[(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)]
+ 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
- )
+ for j in m.set_J[(g,l_star,l)] # only directed arcs
+ ) == 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
- )
-
- # 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
- )
+
+ # 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)
+
+ # *************************************************************************
+ # *************************************************************************
# *****************************************************************************
# *****************************************************************************
\ No newline at end of file
diff --git a/tests/test_esipp_prices.py b/tests/test_esipp_prices.py
index 39dbeb938f72ec5d6dba80d6b35dba247d88a877..813fc0eb311b712ef6b4b85f62869e6e1cc45f74 100644
--- a/tests/test_esipp_prices.py
+++ b/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),