From 093dd504ad1f9fc357577a9d89d60dffcc213dcf Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Pedro=20L=2E=20Magalh=C3=A3es?= <pmlpm@posteo.de>
Date: Wed, 7 Aug 2024 00:28:38 +0200
Subject: [PATCH] Added support for convex and nonconvex price functions using
the lambda formulation, for now only using blocks. Tests were slightly
adjusted too.
---
.../esipp/blocks/{lambda.py => lambda_.py} | 159 ++++++++----------
src/topupopt/problems/esipp/blocks/prices.py | 5 +-
src/topupopt/problems/esipp/utils.py | 4 +-
tests/test_esipp_prices.py | 3 +
4 files changed, 75 insertions(+), 96 deletions(-)
rename src/topupopt/problems/esipp/blocks/{lambda.py => lambda_.py} (75%)
diff --git a/src/topupopt/problems/esipp/blocks/lambda.py b/src/topupopt/problems/esipp/blocks/lambda_.py
similarity index 75%
rename from src/topupopt/problems/esipp/blocks/lambda.py
rename to src/topupopt/problems/esipp/blocks/lambda_.py
index c09943f..7cb6666 100644
--- a/src/topupopt/problems/esipp/blocks/lambda.py
+++ b/src/topupopt/problems/esipp/blocks/lambda_.py
@@ -25,11 +25,6 @@ def price_lambda_block(
# 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)
# *********************************************************************
# *********************************************************************
@@ -53,66 +48,88 @@ def price_lambda_block(
# *********************************************************************
# *********************************************************************
- # 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
- )
+ # set of points
+ b.set_P = 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,
+ within=pyo.NonNegativeReals,
+ initialize=init_param_x_p
+ )
+
+ # 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]
+ # sum(b.param_p_s[i] for i in range(p+1))
+ b.param_y_p = pyo.Param(
+ b.set_P,
+ within=pyo.NonNegativeReals,
+ initialize=init_param_y_p
+ )
- # 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
+ # interpoint weights
+ b.var_weights_p = pyo.Var(b.set_P, within=pyo.UnitInterval)
+
+ # TODO: consider replacing None in the bounds with inf in the parameter
+
+ # transshipment 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_s = pyo.Var(
- b.set_S,
- within=pyo.NonNegativeReals,
- bounds=bounds_var_trans_flows_s
- )
+ b.var_trans_flows = pyo.Var(bounds=bounds_var_trans_flows)
+
+ # *********************************************************************
+ # *********************************************************************
+
+ def rule_constr_stick_to_line(b):
+ return sum(b.var_weights_p[p] for p in b.set_P) == 1
+ b.constr_stick_to_line = pyo.Constraint(rule=rule_constr_stick_to_line)
# *********************************************************************
# *********************************************************************
- # import flow costs and export flow revenues
- def rule_constr_trans_monetary_flows(b):
+ # import flow costs and export flow revenues (y equation)
+ def rule_constr_y_equation(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)
+ sum(b.var_weights_p[p]*b.param_y_p[p] for p in b.set_P)
== 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)
+ sum(b.var_weights_p[p]*b.param_y_p[p] for p in b.set_P)
== b.parent_block().var_efr_glqpk[(g,l,q,p,k)]
)
- b.constr_trans_monetary_flows = pyo.Constraint(
- rule=rule_constr_trans_monetary_flows
- )
+ b.constr_y_equation = pyo.Constraint(rule=rule_constr_y_equation)
+
+ # imported and exported flows (x equation)
+ 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)
+ == b.var_trans_flows
+ )
+ else:
+ return (
+ sum(b.var_weights_p[p]*b.param_x_p[p] for p in b.set_P)
+ == b.var_trans_flows
+ )
+ b.constr_x_equation = pyo.Constraint(rule=rule_constr_x_equation)
- # imported and exported flows
- def rule_constr_trans_flows(b):
+ # imported and exported flows (system equation)
+ def rule_constr_sys_equation(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)
+ ) == b.var_trans_flows
else:
return sum(
b.parent_block().var_v_glljqk[(g,l_star,l,j,q,k)]
@@ -120,8 +137,8 @@ def price_lambda_block(
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)
+ ) == b.var_trans_flows
+ b.constr_sys_equation = pyo.Constraint(rule=rule_constr_sys_equation)
# *********************************************************************
# *********************************************************************
@@ -129,52 +146,10 @@ def price_lambda_block(
# 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
- )
+ # declare SOS2
+ def rule_constr_sos2_weights(b):
+ return [b.var_weights_p[p] for p in b.set_P]
+ b.constr_sos2_weights = pyo.SOSConstraint(rule=rule_constr_sos2_weights, sos=2)
# *********************************************************************
# *********************************************************************
diff --git a/src/topupopt/problems/esipp/blocks/prices.py b/src/topupopt/problems/esipp/blocks/prices.py
index f4f9469..4951fc7 100644
--- a/src/topupopt/problems/esipp/blocks/prices.py
+++ b/src/topupopt/problems/esipp/blocks/prices.py
@@ -3,6 +3,7 @@
import pyomo.environ as pyo
from .other import price_other_block, price_other_no_block
from .delta import price_delta_block, price_delta_no_block
+from .lambda_ import price_lambda_block, price_lambda_no_block
# *****************************************************************************
# *****************************************************************************
@@ -29,7 +30,7 @@ def add_price_functions(
if use_blocks:
# with blocks
if node_price_model == NODE_PRICE_LAMBDA:
- raise NotImplementedError
+ price_lambda_block(model, **kwargs)
elif node_price_model == NODE_PRICE_DELTA:
price_delta_block(model, **kwargs)
else:
@@ -37,7 +38,7 @@ def add_price_functions(
else:
# without blocks
if node_price_model == NODE_PRICE_LAMBDA:
- raise NotImplementedError
+ price_lambda_no_block(model, **kwargs)
elif node_price_model == NODE_PRICE_DELTA:
price_delta_no_block(model, **kwargs)
else:
diff --git a/src/topupopt/problems/esipp/utils.py b/src/topupopt/problems/esipp/utils.py
index bffc0b9..9c5cdbd 100644
--- a/src/topupopt/problems/esipp/utils.py
+++ b/src/topupopt/problems/esipp/utils.py
@@ -98,7 +98,7 @@ def statistics(ipp: InfrastructurePlanningProblem,
)
if prices_in_block:
# imports
- if ipp.node_price_model == NODE_PRICE_DELTA:
+ if ipp.node_price_model == NODE_PRICE_DELTA or ipp.node_price_model == NODE_PRICE_LAMBDA:
# imports
imports_qpk = {
qpk: pyo.value(
@@ -152,7 +152,7 @@ def statistics(ipp: InfrastructurePlanningProblem,
}
else:
# not in a block
- if ipp.node_price_model == NODE_PRICE_DELTA:
+ if ipp.node_price_model == NODE_PRICE_DELTA or ipp.node_price_model == NODE_PRICE_LAMBDA:
# imports
imports_qpk = {
qpk: pyo.value(
diff --git a/tests/test_esipp_prices.py b/tests/test_esipp_prices.py
index 2652bd6..39dbeb9 100644
--- a/tests/test_esipp_prices.py
+++ b/tests/test_esipp_prices.py
@@ -221,6 +221,7 @@ class TestESIPPProblem:
# no sos, regular time intervals
ipp = build_solve_ipp(
+ solver='scip' if node_price_model == NODE_PRICE_LAMBDA else 'glpk',
solver_options={},
perform_analysis=False,
plot_results=False, # True,
@@ -351,6 +352,7 @@ class TestESIPPProblem:
# no sos, regular time intervals
ipp = build_solve_ipp(
+ solver='scip' if node_price_model == NODE_PRICE_LAMBDA else 'glpk',
solver_options={},
perform_analysis=False,
plot_results=False, # True,
@@ -680,6 +682,7 @@ class TestESIPPProblem:
# no sos, regular time intervals
ipp = build_solve_ipp(
+ solver='scip' if node_price_model == NODE_PRICE_LAMBDA else 'glpk',
solver_options={},
perform_analysis=False,
plot_results=False, # True,
--
GitLab