Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • pmag/topupopt
1 result
Show changes
Commits on Source (18)
Showing
with 3248 additions and 2290 deletions
# -*- coding: utf-8 -*-
# from . import mvesipp
# -*- coding: utf-8 -*-
\ No newline at end of file
# -*- coding: utf-8 -*-
\ No newline at end of file
# 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()
# *********************************************************************
# *********************************************************************
# 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
)
# segment usage variables
b.var_segment_usage_s = pyo.Var(
b.set_S,
within=pyo.UnitInterval,
)
# *********************************************************************
# *********************************************************************
# import flow costs and export flow revenues (y function)
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]
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]
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 (x function)
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 defines the x value
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
# variables to indicate if the segments are active
b.var_segment_is_full = pyo.Var(
b.set_S,
within=pyo.Binary
)
# if next level is active, bin. var. for current level must be one
# if bin. var. for current level is one, segment must be used fully
# if previous level is not full, bin. var for current level must be zero
# if bin. var. for current level is zero,
def rule_constr_nonconvex_p1(b, s):
if s == len(b.set_S)-1:
# last segment, skip
return pyo.Constraint.Skip
return (
b.var_segment_usage_s[s+1] <=
b.var_segment_is_full[s]
)
b.constr_nonconvex_p1 = pyo.Constraint(
b.set_S, rule=rule_constr_nonconvex_p1
)
def rule_constr_nonconvex_p2(b, s):
if s == len(b.set_S)-1:
# last segment, skip
return pyo.Constraint.Skip
return (
b.var_segment_usage_s[s] >=
b.var_segment_is_full[s]
)
b.constr_nonconvex_p2 = pyo.Constraint(
b.set_S, rule=rule_constr_nonconvex_p2
)
def rule_constr_nonconvex_p3(b, s):
if s == len(b.set_S)-1:
# last segment, skip
return pyo.Constraint.Skip
return (
b.var_segment_usage_s[s] >=
b.var_segment_is_full[s+1]
)
b.constr_nonconvex_p3 = pyo.Constraint(
b.set_S, rule=rule_constr_nonconvex_p3
)
def rule_constr_nonconvex_p4(b, s):
if s == len(b.set_S)-1:
# last segment, skip
return pyo.Constraint.Skip
return (
b.var_segment_usage_s[s+1] <=
b.var_segment_is_full[s+1]
)
b.constr_nonconvex_p4 = pyo.Constraint(
b.set_S, rule=rule_constr_nonconvex_p4
)
# *********************************************************************
# *********************************************************************
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
):
# 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_glqpk(m, g, l, q, p, k):
try:
return (0, sum(m.param_v_max_glqpks[s] for s in m.set_S[(g, l, q, p, k)]))
except Exception:
return (0, None)
model.var_trans_flows_glqpk = pyo.Var(
model.set_GLQPK, within=pyo.NonNegativeReals, bounds=bounds_var_trans_flows_glqpk
)
# segment usage variables
model.var_segment_usage_glqpks = pyo.Var(
model.set_GLQPKS,
within=pyo.UnitInterval,
)
# *************************************************************************
# *************************************************************************
# import flow costs and export flow revenues (y function)
def rule_constr_trans_monetary_flows(m, g, l, q, p, k):
if (g,l) in m.set_GL_imp:
return (
sum(
m.var_segment_usage_glqpks[(g, l, q, p, k, s)]*
m.param_p_glqpks[(g, l, q, p, k, s)]*
m.param_v_max_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_segment_usage_glqpks[(g, l, q, p, k, s)]*
m.param_p_glqpks[(g, l, q, p, k, s)]*
m.param_v_max_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 (x function)
def rule_constr_trans_flows_seg(m, g, l, q, p, k):
return sum(
m.var_segment_usage_glqpks[(g, l, q, p, k, s)]*m.param_v_max_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)]
model.constr_trans_flows_seg = pyo.Constraint(
model.set_GLQPK,
rule=rule_constr_trans_flows_seg
)
# imported and exported flows: energy system defines the x value
def rule_constr_trans_flows_sys(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
) == 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)]
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
) == m.var_trans_flows_glqpk[(g, l, q, p, k)]
model.constr_trans_flows_sys = pyo.Constraint(
model.set_GLQPK,
rule=rule_constr_trans_flows_sys
)
# *************************************************************************
# *************************************************************************
# non-convex price functions
# variables to indicate if the segments are active
model.var_segment_is_full = pyo.Var(
model.set_GLQPKS,
within=pyo.Binary
)
# if next level is active, bin. var. for current level must be one
# if bin. var. for current level is one, segment must be used fully
# if previous level is not full, bin. var for current level must be zero
# if bin. var. for current level is zero,
def rule_constr_nonconvex_p1(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_segment_usage_glqpks[(g,l,q,p,k,s+1)] <=
m.var_segment_is_full[(g,l,q,p,k,s)]
)
model.constr_nonconvex_p1 = pyo.Constraint(
model.set_GLQPKS, rule=rule_constr_nonconvex_p1
)
def rule_constr_nonconvex_p2(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_segment_usage_glqpks[(g,l,q,p,k,s)] >=
m.var_segment_is_full[(g,l,q,p,k,s)]
)
model.constr_nonconvex_p2 = pyo.Constraint(
model.set_GLQPKS, rule=rule_constr_nonconvex_p2
)
def rule_constr_nonconvex_p3(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_segment_usage_glqpks[(g,l,q,p,k,s)] >=
m.var_segment_is_full[(g,l,q,p,k,s+1)]
)
model.constr_nonconvex_p3 = pyo.Constraint(
model.set_GLQPKS, rule=rule_constr_nonconvex_p3
)
def rule_constr_nonconvex_p4(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_segment_usage_glqpks[(g,l,q,p,k,s+1)] <=
m.var_segment_is_full[(g,l,q,p,k,s+1)]
)
model.constr_nonconvex_p4 = pyo.Constraint(
model.set_GLQPKS, rule=rule_constr_nonconvex_p4
)
# *****************************************************************************
# *****************************************************************************
\ No newline at end of file
# 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()
# *********************************************************************
# *********************************************************************
# 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
# *********************************************************************
# *********************************************************************
# set of points
b.set_O = pyo.Set(initialize=[i for i in range(len(b.set_S)+1)])
# set of x points
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_o
)
# set of y points
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_o = pyo.Param(
b.set_O,
within=pyo.NonNegativeReals,
initialize=init_param_y_o
)
# interpoint weights
b.var_weights_o = pyo.Var(b.set_O, 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 = pyo.Var(bounds=bounds_var_trans_flows)
# *********************************************************************
# *********************************************************************
def rule_constr_stick_to_line(b):
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)
# *********************************************************************
# *********************************************************************
# 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_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_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)
# 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_o[p]*b.param_x_o[p] for p in b.set_O)
== b.var_trans_flows
)
else:
return (
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)
# 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
) == 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_sys_equation = pyo.Constraint(rule=rule_constr_sys_equation)
# *********************************************************************
# *********************************************************************
# non-convex price functions
if not b.param_price_function_is_convex:
# declare SOS2
def rule_constr_sos2_weights(b):
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)
# *********************************************************************
# *********************************************************************
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
):
# 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
# *************************************************************************
# *************************************************************************
# 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_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 (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_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_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_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 (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
) == 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)]
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
) == 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
# 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
# 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
......@@ -3,7 +3,7 @@
import pyomo.environ as pyo
from .blocks.networks import add_network_restrictions
from .blocks.prices import add_prices_block
from .blocks.prices import add_price_functions
# *****************************************************************************
# *****************************************************************************
......@@ -11,9 +11,11 @@ from .blocks.prices import add_prices_block
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,
enable_initialisation: bool = True
):
# TODO: test initialisation
......@@ -57,33 +59,26 @@ def create_model(
# *************************************************************************
# set of assessments
model.set_Q = pyo.Set()
# TODO: use rangesets for time-related sets
# set of time step intervals for each assessment
model.set_K_q = pyo.Set(model.set_Q)
# set of representative evaluation periods for each assessment
model.set_P_q = pyo.Set(model.set_Q)
# set of networks
model.set_G = pyo.Set()
# set of nodes on each network
model.set_L = pyo.Set(model.set_G)
# set of importing nodes on each network
model.set_L_imp = pyo.Set(model.set_G, within=model.set_L)
# set of exporting nodes on each network
model.set_L_exp = pyo.Set(model.set_G, within=model.set_L)
# *************************************************************************
......@@ -963,20 +958,6 @@ def create_model(
initialize=(init_set_GLLJH_arc_inv_sos1 if enable_initialisation else None),
)
# set of GLLJ-indexed GLLJH tuples for new arcs modelled using SOS1
def init_set_GLLJH_arc_inv_sos1_gllj(m, g, l1, l2, j):
return ((g, l1, l2, j, h) for h in m.set_H_gllj[(g, l1, l2, j)])
model.set_GLLJH_arc_inv_sos1_gllj = pyo.Set(
model.set_GLLJ_arc_inv_sos1,
dimen=5,
within=model.set_GLLJH_arc_inv_sos1,
initialize=(
init_set_GLLJH_arc_inv_sos1_gllj if enable_initialisation else None
),
)
# *************************************************************************
# flow sense determination using SOS1
......@@ -1315,18 +1296,6 @@ def create_model(
dimen=2, within=model.set_TH, initialize=init_set_TH_arc_inv_sos1
)
# set of t-indexed TH tuples for groups of arcs relying on SOS1
def init_set_TH_arc_inv_sos1_t(m, t):
return ((t, h) for h in m.set_H_t[t])
model.set_TH_arc_inv_sos1_t = pyo.Set(
model.set_T_sos1,
dimen=2,
within=model.set_TH_arc_inv_sos1,
initialize=init_set_TH_arc_inv_sos1_t,
)
# minimum cost of a group of arcs
model.param_c_arc_min_th = pyo.Param(model.set_TH, within=pyo.NonNegativeReals)
......@@ -2163,7 +2132,11 @@ def create_model(
)
# prices
add_prices_block(model)
add_price_functions(
model,
use_blocks=use_prices_block,
node_price_model=node_price_model
)
# *************************************************************************
# *************************************************************************
......@@ -2435,14 +2408,21 @@ def create_model(
# *************************************************************************
# SOS1 constraints for arc selection
# SOS1 constraints for arc selection
def rule_constr_arc_sos1(m, g, l1, l2, j):
var_list = [
m.var_delta_arc_inv_glljh[(g, l1, l2, j, h)]
for h in m.set_H_gllj[(g, l1, l2, j)]
]
weight_list = [
m.param_arc_inv_sos1_weights_glljh[(g, l1, l2, j, h)]
for h in m.set_H_gllj[(g, l1, l2, j)]
]
return (var_list, weight_list)
model.constr_arc_sos1 = pyo.SOSConstraint(
model.set_GLLJ_arc_inv_sos1, # for (directed and undirected) new arcs
var=model.var_delta_arc_inv_glljh, # set_GLLJH_sgl indexes the variables
index=model.set_GLLJH_arc_inv_sos1_gllj, # key: GLLJ; value: GLLJH
weights=model.param_arc_inv_sos1_weights_glljh, # key: GLLJH; alue: weight
sos=1,
model.set_GLLJ_arc_inv_sos1,
rule=rule_constr_arc_sos1,
sos=1
)
# *************************************************************************
......@@ -2572,16 +2552,29 @@ def create_model(
# *************************************************************************
# SOS1 constraints for flow sense determination (undirected arcs)
# model.constr_sns_sos1 = pyo.SOSConstraint(
# model.set_GLLJ_und, # one constraint per undirected arc
# model.set_QK, # and time interval
# var=model.var_zeta_sns_glljqk, # set_GLLJ_und_red and set_K
# index=model.set_GLLJQK_und_sns_sos1_red_gllj, #
# weights=model.param_arc_sns_sos1_weights_glljqk,
# sos=1,
# )
# TODO: make the weight list dependent in model options rather than literal
def rule_constr_sns_sos1(m, g, l1, l2, j, q, k):
var_list = [
m.var_zeta_sns_glljqk[(g, l1, l2, j, q, k)],
m.var_zeta_sns_glljqk[(g, l2, l1, j, q, k)]
]
weight_list = [1, 2] if True else [2, 1]
return (var_list, weight_list)
model.constr_sns_sos1 = pyo.SOSConstraint(
model.set_GLLJ_und, # one constraint per undirected arc
model.set_QK, # and time interval
var=model.var_zeta_sns_glljqk, # set_GLLJ_und_red and set_K
index=model.set_GLLJQK_und_sns_sos1_red_gllj, #
weights=model.param_arc_sns_sos1_weights_glljqk,
sos=1,
model.set_GLLJ_und,
model.set_QK,
rule=rule_constr_sns_sos1,
sos=1
)
# *************************************************************************
# static losses
......@@ -2941,13 +2934,20 @@ def create_model(
# *************************************************************************
# SOS1 constraints for arc group selection
def rule_constr_arc_group_sos1(m, t):
var_list = [
m.var_delta_arc_inv_th[(t, h)]
for h in m.set_H_t[t]
]
weight_list = [
m.param_arc_inv_sos1_weights_th[(t, h)]
for h in m.set_H_t[t]
]
return (var_list, weight_list)
model.constr_arc_group_sos1 = pyo.SOSConstraint(
model.set_T_sos1, # for all groups using sos1
var=model.var_delta_arc_inv_th, # set_TH indexes the variables
index=model.set_TH_arc_inv_sos1_t, # key: t; value: TH
weights=model.param_arc_inv_sos1_weights_th, # key: TH; value: weight
sos=1,
model.set_T_sos1,
rule=rule_constr_arc_group_sos1,
sos=1
)
# *************************************************************************
......
This diff is collapsed.
......@@ -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,8 @@ 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
# *********************************************************************
......@@ -122,6 +115,16 @@ class InfrastructurePlanningProblem(EnergySystem):
# self.number_time_intervals = {
# q: len(time_frame.time_intervals[q]) for q in self.time_frame.assessments
# }
# *********************************************************************
# *********************************************************************
# options
self.use_prices_block = use_prices_block
self.node_price_model = node_price_model
# *********************************************************************
# *********************************************************************
# initialise
......@@ -470,6 +473,7 @@ class InfrastructurePlanningProblem(EnergySystem):
use_sos1: bool = False,
use_interface: bool = False,
use_nnr_variables_if_possible: bool = False,
sos1_weight_method: int = SOS1_ARC_WEIGHTS_NONE
) -> int:
"""
Create a group of arcs whose invesment is to be decided collectively.
......@@ -490,14 +494,14 @@ class InfrastructurePlanningProblem(EnergySystem):
Returns the key that identifies the arc group, an integer.
"""
# make sure there is at least one arc
if len(gllj_tuples) < 2:
raise ValueError(
"At least two arcs need to be identified to create a group."
)
for arc_number, gllj in enumerate(gllj_tuples):
# does the network exist?
......@@ -555,10 +559,10 @@ class InfrastructurePlanningProblem(EnergySystem):
# groups using interfaces
self.groups_int[new_t] = use_interface
# TODO: allow users to set the weights
# groups using sos1 for arc selection
self.groups_arc_sos1[new_t] = use_sos1
if use_sos1:
self.groups_arc_sos1[new_t] = [i for i in range(number_options)]
# groups using nnr for arc selection
......@@ -1625,28 +1629,22 @@ class InfrastructurePlanningProblem(EnergySystem):
raise ValueError(
"The method to determine the SOS1 weights was not recognised."
)
if method == self.SOS1_ARC_WEIGHTS_CAP:
sos1_weights = tuple(capacity for capacity in arcs.capacity)
elif method == self.SOS1_ARC_WEIGHTS_COST:
sos1_weights = tuple(cost for cost in arcs.minimum_cost)
elif method == self.SOS1_ARC_WEIGHTS_SPEC_CAP:
sos1_weights = tuple(
cap / cost for cap, cost in zip(arcs.capacity, arcs.minimum_cost)
)
elif method == self.SOS1_ARC_WEIGHTS_SPEC_COST:
sos1_weights = tuple(
cost / cap for cap, cost in zip(arcs.capacity, arcs.minimum_cost)
)
else: # SOS1_ARC_WEIGHTS_NONE
return None
# make sure they are unique
if verify_weights:
for weight in sos1_weights:
if sos1_weights.count(weight) >= 2: # TODO: reach this point
......@@ -1655,7 +1653,6 @@ class InfrastructurePlanningProblem(EnergySystem):
+ "special ordered sets of type 1 (SOS1),"
+ " since some weights are not unique."
)
return sos1_weights
# *************************************************************************
......@@ -1663,10 +1660,12 @@ 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)
self.model = create_model(
name,
use_prices_block=self.use_prices_block,
node_price_model=self.node_price_model
)
# *************************************************************************
# *************************************************************************
......@@ -1883,19 +1882,6 @@ class InfrastructurePlanningProblem(EnergySystem):
.number_segments()
)
)
if not self.networks[g].nodes[l][Network.KEY_NODE_PRICES_TIME_INVARIANT]
else tuple(
s
for s in range(
self.networks[g]
.nodes[l][Network.KEY_NODE_PRICES][(q, p, k)]
.number_segments()
)
)
# for g in self.networks.keys()
# for l in self.networks[g].nodes
# if (l in self.networks[g].import_nodes or
# l in self.networks[g].export_nodes)
for (g, l) in set_GL_exp_imp
for (q, p, k) in set_QPK
}
......@@ -2467,7 +2453,7 @@ class InfrastructurePlanningProblem(EnergySystem):
# set of arc groups relying on SOS1
set_T_sos1 = tuple(t for t in set_T if self.groups_arc_sos1[t])
set_T_sos1 = tuple(t for t in set_T if t in self.groups_arc_sos1)
# set of arg groups relying on non-negative real variables
......@@ -2561,7 +2547,7 @@ class InfrastructurePlanningProblem(EnergySystem):
}
# maximum resource volume per segment (infinity is the default)
param_v_max_glqpks = {
(g, l, q, p, k, s): self.networks[g]
.nodes[l][Network.KEY_NODE_PRICES][(q, p, k)]
......@@ -3145,12 +3131,6 @@ class InfrastructurePlanningProblem(EnergySystem):
set_GLLJH_arc_inv_sos1 = tuple(param_arc_inv_sos1_weights_glljh.keys())
set_GLLJH_arc_inv_sos1_gllj = {
(g, u, v, j): tuple((g, u, v, j, h) for h in set_H_gllj[(g, u, v, j)])
for (g, u, v) in set_J_arc_sos1
for j in set_J_arc_sos1[(g, u, v)]
}
# sos1 weights for flow sense determination
param_arc_sns_sos1_weights_glljqk = {
......@@ -3214,7 +3194,7 @@ class InfrastructurePlanningProblem(EnergySystem):
}
param_arc_inv_sos1_weights_th = {
(t, h): h
(t, h): self.groups_arc_sos1[t][h]
# (t, h): mean(
# (self.use_sos1_arc_inv[(g,u,v,j)][h]
# if self.use_sos1_arc_inv[(g,u,v,j)] is not None else h)
......@@ -3309,6 +3289,8 @@ class InfrastructurePlanningProblem(EnergySystem):
# *********************************************************************
# produce a dictionary with the data for the problem
# TODO: built the dict as a function of the price submodel
data_dict = {
None: {
......@@ -3433,7 +3415,6 @@ class InfrastructurePlanningProblem(EnergySystem):
"set_GLLJ_pre_fin_red": set_GLLJ_pre_fin_red,
"set_GLLJ_arc_inv_sos1": set_GLLJ_arc_inv_sos1,
"set_GLLJH_arc_inv_sos1": set_GLLJH_arc_inv_sos1,
"set_GLLJH_arc_inv_sos1_gllj": set_GLLJH_arc_inv_sos1_gllj,
"set_GLLJQK_und_sns_sos1_red": set_GLLJQK_und_sns_sos1_red,
"set_GLLJ_int": set_GLLJ_int,
"set_GLLJ_static": set_GLLJ_static,
......@@ -3524,6 +3505,16 @@ class InfrastructurePlanningProblem(EnergySystem):
"param_arc_sns_sos1_weights_glljqk": param_arc_sns_sos1_weights_glljqk,
# *****************************************************************
# *****************************************************************
"block_prices": {
(*gl,*qpk): {
'param_price_function_is_convex': {None: param_price_function_is_convex[(*gl,*qpk)]},
'set_S': {None: set_S[(*gl,*qpk)]},
'param_p_s': {s: param_p_glqpks[(*gl,*qpk,s)] for s in set_S[(*gl,*qpk)]},
'param_v_max_s': {s: param_v_max_glqpks[(*gl,*qpk,s)] for s in set_S[(*gl,*qpk)] if (*gl,*qpk,s) in param_v_max_glqpks},
}
for gl in set_GL_exp_imp
for qpk in set_QPK
}
}
}
......
......@@ -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
# *****************************************************************************
# *****************************************************************************
......@@ -69,6 +70,7 @@ def statistics(ipp: InfrastructurePlanningProblem,
other_node_keys: tuple = None):
"Returns flow statistics using the optimisation results."
prices_in_block = ipp.use_prices_block
if type(import_node_keys) == type(None):
# pick all import nodes
......@@ -94,36 +96,116 @@ def statistics(ipp: InfrastructurePlanningProblem,
for node_key in net.nodes()
if Network.KEY_NODE_BASE_FLOW in net.nodes[node_key]
)
# imports
imports_qpk = {
qpk: pyo.value(
sum(
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
for s in ipp.instance.set_S[(g,l_imp,*qpk)]
)
*ipp.instance.param_c_time_qpk[qpk]
)
for qpk in ipp.time_frame.qpk()
}
# exports
exports_qpk = {
qpk: pyo.value(
sum(
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
for s in ipp.instance.set_S[(g,l_exp,*qpk)]
)
*ipp.instance.param_c_time_qpk[qpk]
)
for qpk in ipp.time_frame.qpk()
}
if prices_in_block:
# imports
if ipp.node_price_model == NODE_PRICE_DELTA or ipp.node_price_model == NODE_PRICE_LAMBDA:
# 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]
)
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]
)
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
if ipp.node_price_model == NODE_PRICE_DELTA or ipp.node_price_model == NODE_PRICE_LAMBDA:
# imports
imports_qpk = {
qpk: pyo.value(
sum(
ipp.instance.var_trans_flows_glqpk[(g,l_imp,*qpk)]
for g, l_imp in import_node_keys
)
*ipp.instance.param_c_time_qpk[qpk]
)
for qpk in ipp.time_frame.qpk()
}
# exports
exports_qpk = {
qpk: pyo.value(
sum(
ipp.instance.var_trans_flows_glqpk[(g,l_exp,*qpk)]
for g, l_exp in export_node_keys
)
*ipp.instance.param_c_time_qpk[qpk]
)
for qpk in ipp.time_frame.qpk()
}
else:
# imports
imports_qpk = {
qpk: pyo.value(
sum(
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
for s in ipp.instance.set_S[(g,l_imp,*qpk)]
)
*ipp.instance.param_c_time_qpk[qpk]
)
for qpk in ipp.time_frame.qpk()
}
# exports
exports_qpk = {
qpk: pyo.value(
sum(
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
for s in ipp.instance.set_S[(g,l_exp,*qpk)]
)
*ipp.instance.param_c_time_qpk[qpk]
)
for qpk in ipp.time_frame.qpk()
}
# balance
balance_qpk = {
qpk: imports_qpk[qpk]-exports_qpk[qpk]
......
This diff is collapsed.
......@@ -651,7 +651,7 @@ class TestDataFinance:
# trigger ValueError
error_triggered = False
error_raised = False
investment_period = analysis_period_span + 1
try:
npv_salvage = present_salvage_value_annuity(
......@@ -662,8 +662,8 @@ class TestDataFinance:
analysis_period_span=analysis_period_span,
)
except ValueError:
error_triggered = True
assert error_triggered
error_raised = True
assert error_raised
# *********************************************************************
......@@ -1197,7 +1197,7 @@ class TestDataFinance:
investment_period = analysis_period_span + 1
error_triggered = False
error_raised = False
try:
residual_value = salvage_value_linear_depreciation(
investment=investment,
......@@ -1206,8 +1206,8 @@ class TestDataFinance:
analysis_period_span=analysis_period_span,
)
except ValueError:
error_triggered = True
assert error_triggered
error_raised = True
assert error_raised
# *************************************************************************
# *************************************************************************
......@@ -1318,78 +1318,78 @@ class TestDataFinance:
# TypeError('The discount rates must be provided as a tuple.')
error_triggered = False
error_raised = False
try:
my_inv = Investment(list(i_t), R_t)
except TypeError:
error_triggered = True
assert error_triggered
error_raised = True
assert error_raised
# *********************************************************************
# ValueError('The duration of the period under analysis must be positive.')
error_triggered = False
error_raised = False
try:
my_inv = Investment(tuple())
except ValueError:
error_triggered = True
assert error_triggered
error_raised = True
assert error_raised
# *********************************************************************
# TypeError('The discount rate must be provided as a float.')
error_triggered = False
error_raised = False
try:
my_inv = Investment(None, None, 5, 10)
except TypeError:
error_triggered = True
assert error_triggered
error_raised = True
assert error_raised
# *********************************************************************
# ValueError('The discount rate must be in the open interval between 0 and 1.)
error_triggered = False
error_raised = False
try:
my_inv = Investment(None, None, 1.35, 10)
except ValueError:
error_triggered = True
assert error_triggered
error_raised = True
assert error_raised
# *********************************************************************
# TypeError('The duration of the period under consideration must be provided as an integer.')
error_triggered = False
error_raised = False
try:
my_inv = Investment(None, None, 0.35, 10.0)
except TypeError:
error_triggered = True
assert error_triggered
error_raised = True
assert error_raised
# *********************************************************************
# ValueError('The duration of the period under analysis must be positive.)
error_triggered = False
error_raised = False
try:
my_inv = Investment(None, None, 0.35, 0)
except ValueError:
error_triggered = True
assert error_triggered
error_raised = True
assert error_raised
# *********************************************************************
# TypeError('The net cash flows must be provided as a list.')
error_triggered = False
error_raised = False
try:
my_inv = Investment(i_t, tuple(R_t))
except TypeError:
error_triggered = True
assert error_triggered
error_raised = True
assert error_raised
# *********************************************************************
......
......@@ -140,7 +140,7 @@ class TestDataUtils:
# raise exception
error_triggered = False
error_raised = False
time_interval_durations.pop(0)
try:
new_profile = utils.create_profile_using_time_weighted_state(
......@@ -151,8 +151,8 @@ class TestDataUtils:
states_correlate_profile=True,
)
except ValueError:
error_triggered = True
assert error_triggered
error_raised = True
assert error_raised
# *********************************************************************
# *********************************************************************
......@@ -333,14 +333,14 @@ class TestDataUtils:
# use zero iterations to force an error
error_triggered = False
error_raised = False
try:
new_key = utils.generate_pseudo_unique_key(
key_list=key_list, max_iterations=0
)
except Exception:
error_triggered = True
assert error_triggered
error_raised = True
assert error_raised
# use a seed number to trigger more iterations
......
......@@ -654,7 +654,7 @@ class TestDistrictHeatingNetwork:
# # single pipe, external cost, offset
# error_triggered = False
# error_raised = False
# try:
# pipe_trench_obj = PipeTrench(name='hello',
# trenches={0: trench_tech},
......@@ -665,12 +665,12 @@ class TestDistrictHeatingNetwork:
# minimum_cost_offset=external_cost,
# validate=True)
# except TypeError:
# error_triggered = True
# assert error_triggered
# error_raised = True
# assert error_raised
# # use list as minimum cost offset
# error_triggered = False
# error_raised = False
# try:
# pipe_trench_obj = PipeTrench(name='hello',
# trenches={0: trench_tech},
......@@ -684,8 +684,8 @@ class TestDistrictHeatingNetwork:
# ),
# validate=True)
# except TypeError:
# error_triggered = True
# assert error_triggered
# error_raised = True
# assert error_raised
# #**************************************************************************
# #**************************************************************************
......@@ -754,7 +754,7 @@ class TestDistrictHeatingNetwork:
# # single pipe, external cost, offset
# error_triggered = False
# error_raised = False
# try:
# pipe_trench_obj = PipeTrench(name='hello',
# trenches={0: trench_tech},
......@@ -765,12 +765,12 @@ class TestDistrictHeatingNetwork:
# minimum_cost_offset=external_cost,
# validate=True)
# except TypeError:
# error_triggered = True
# assert error_triggered
# error_raised = True
# assert error_raised
# # use list as minimum cost offset
# error_triggered = False
# error_raised = False
# try:
# pipe_trench_obj = PipeTrench(name='hello',
# trenches={0: trench_tech},
......@@ -784,8 +784,8 @@ class TestDistrictHeatingNetwork:
# ),
# validate=True)
# except TypeError:
# error_triggered = True
# assert error_triggered
# error_raised = True
# assert error_raised
# #**************************************************************************
# #**************************************************************************
......
......@@ -20,7 +20,6 @@ from topupheat.common.fluids import FluidDatabase # , Fluid
# *****************************************************************************
# *****************************************************************************
class TestDistrictHeatingNetworkUtils:
# *************************************************************************
# *************************************************************************
......@@ -528,7 +527,7 @@ class TestDistrictHeatingNetworkUtils:
# update the nodes
network.add_node(0, x=55, y=12)
network.add_node(2, x=55.01, y=12.01)
# *********************************************************************
utils.summarise_network_by_pipe_technology(network, False)
......
# imports
# local
# import numpy as np
# import networkx as nx
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
# *****************************************************************************
# *****************************************************************************
def check_problem_size(ipp: InfrastructurePlanningProblem, nc, nv, nnz):
assert ipp.results["Problem"][0]["Number of constraints"] == nc # should be 80
assert ipp.results["Problem"][0]["Number of variables"] == nv # should be 84
assert ipp.results["Problem"][0]["Number of nonzeros"] == nnz
# *****************************************************************************
# *****************************************************************************
def build_solve_ipp(
solver: str = 'glpk',
solver_options: dict = None,
use_sos_arcs: bool = False,
use_sos_arc_groups: bool = False,
arc_sos_weight_key: str = InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_NONE,
arc_use_real_variables_if_possible: bool = False,
use_sos_sense: bool = False,
sense_sos_weight_key: int = (
InfrastructurePlanningProblem.SOS1_SENSE_WEIGHT_NOMINAL_HIGHER
),
sense_use_real_variables_if_possible: bool = False,
sense_use_arc_interfaces: bool = False,
perform_analysis: bool = False,
plot_results: bool = False,
print_solver_output: bool = False,
time_frame: EconomicTimeFrame = None,
networks: dict = None,
converters: dict = None,
static_losses_mode=None,
mandatory_arcs: list = None,
max_number_parallel_arcs: dict = None,
arc_groups_dict: dict = None,
init_aux_sets: bool = False,
# discount_rates: dict = None,
assessment_weights: dict = None,
simplify_problem: bool = False,
use_prices_block: bool = False,
node_price_model: int = NODE_PRICE_OTHER
):
if type(assessment_weights) != dict:
assessment_weights = {} # default
if type(converters) != dict:
converters = {}
# time weights
# relative weight of time period
# one interval twice as long as the average is worth twice
# one interval half as long as the average is worth half
# time_weights = [
# [time_period_duration/average_time_interval_duration
# for time_period_duration in intraperiod_time_interval_duration]
# for p in range(number_periods)]
time_weights = None # nothing yet
normalised_time_interval_duration = None # nothing yet
# create problem object
ipp = InfrastructurePlanningProblem(
# discount_rates=discount_rates,
time_frame=time_frame,
# reporting_periods=time_frame.reporting_periods,
# time_intervals=time_frame.time_interval_durations,
time_weights=time_weights,
normalised_time_interval_duration=normalised_time_interval_duration,
assessment_weights=assessment_weights,
use_prices_block=use_prices_block,
node_price_model=node_price_model
)
# add networks and systems
for netkey, net in networks.items():
ipp.add_network(network_key=netkey, network=net)
# add converters
for cvtkey, cvt in converters.items():
ipp.add_converter(converter_key=cvtkey, converter=cvt)
# define arcs as mandatory
if type(mandatory_arcs) == list:
for full_arc_key in mandatory_arcs:
ipp.make_arc_mandatory(full_arc_key[0], full_arc_key[1:])
# if make_all_arcs_mandatory:
# for network_key in ipp.networks:
# for arc_key in ipp.networks[network_key].edges(keys=True):
# # preexisting arcs are no good
# if ipp.networks[network_key].edges[arc_key][
# Network.KEY_ARC_TECH].has_been_selected():
# continue
# ipp.make_arc_mandatory(network_key, arc_key)
# set up the use of sos for arc selection
if use_sos_arcs:
for network_key in ipp.networks:
for arc_key in ipp.networks[network_key].edges(keys=True):
if (
ipp.networks[network_key]
.edges[arc_key][Network.KEY_ARC_TECH]
.has_been_selected()
):
# skip arcs that have already been selected (pre-existing)
continue
ipp.use_sos1_for_arc_selection(
network_key,
arc_key,
use_real_variables_if_possible=(
arc_use_real_variables_if_possible
),
sos1_weight_method=arc_sos_weight_key,
)
# set up the use of sos for flow sense determination
if use_sos_sense:
for network_key in ipp.networks:
for arc_key in ipp.networks[network_key].edges(keys=True):
if not ipp.networks[network_key].edges[arc_key][
Network.KEY_ARC_UND
]:
continue
ipp.use_sos1_for_flow_senses(
network_key,
arc_key,
use_real_variables_if_possible=(
sense_use_real_variables_if_possible
),
use_interface_variables=sense_use_arc_interfaces,
sos1_weight_method=sense_sos_weight_key,
)
elif sense_use_arc_interfaces: # set up the use of arc interfaces w/o sos1
for network_key in ipp.networks:
for arc_key in ipp.networks[network_key].edges(keys=True):
if (
ipp.networks[network_key]
.edges[arc_key][Network.KEY_ARC_TECH]
.has_been_selected()
):
continue
ipp.use_interface_variables_for_arc_selection(network_key, arc_key)
# static losses
if static_losses_mode == ipp.STATIC_LOSS_MODE_ARR:
ipp.place_static_losses_arrival_node()
elif static_losses_mode == ipp.STATIC_LOSS_MODE_DEP:
ipp.place_static_losses_departure_node()
elif static_losses_mode == ipp.STATIC_LOSS_MODE_US:
ipp.place_static_losses_upstream()
elif static_losses_mode == ipp.STATIC_LOSS_MODE_DS:
ipp.place_static_losses_downstream()
else:
raise ValueError("Unknown static loss modelling mode.")
# *********************************************************************
# groups
if type(arc_groups_dict) != type(None):
for key in arc_groups_dict:
ipp.create_arc_group(
arc_groups_dict[key],
use_sos1=use_sos_arc_groups,
sos1_weight_method=arc_sos_weight_key
)
# *********************************************************************
# maximum number of parallel arcs
for key in max_number_parallel_arcs:
ipp.set_maximum_number_parallel_arcs(
network_key=key[0],
node_a=key[1],
node_b=key[2],
limit=max_number_parallel_arcs[key],
)
# *********************************************************************
if simplify_problem:
ipp.simplify_peak_total_assessments()
# *********************************************************************
# instantiate (disable the default case v-a-v fixed losses)
# ipp.instantiate(place_fixed_losses_upstream_if_possible=False)
ipp.instantiate(initialise_ancillary_sets=init_aux_sets)
# optimise
ipp.optimise(
solver_name=solver,
solver_options=solver_options,
output_options={},
print_solver_output=print_solver_output,
)
# ipp.instance.pprint()
# return the problem object
return ipp
# *****************************************************************************
# *****************************************************************************
\ No newline at end of file
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
......@@ -598,145 +598,145 @@ class TestResourcePrice:
# create object without prices
error_triggered = False
error_raised = False
try:
_ = ResourcePrice(prices=None, volumes=volumes)
except TypeError:
error_triggered = True
assert error_triggered
error_raised = True
assert error_raised
# *********************************************************************
# create object with negative prices in lists
error_triggered = False
error_raised = False
try:
_ = ResourcePrice(prices=[7, -3, 2], volumes=[3, 4, 5])
except ValueError:
error_triggered = True
assert error_triggered
error_raised = True
assert error_raised
# *********************************************************************
# create object where an intermediate segment has no volume limit
error_triggered = False
error_raised = False
try:
_ = ResourcePrice(prices=[7, 4, 2], volumes=[3, None, 5])
except ValueError:
error_triggered = True
assert error_triggered
error_raised = True
assert error_raised
# *********************************************************************
# create object with negative volumes in lists
error_triggered = False
error_raised = False
try:
_ = ResourcePrice(prices=[7, 3, 2], volumes=[4, -1, 2])
except ValueError:
error_triggered = True
assert error_triggered
error_raised = True
assert error_raised
# *********************************************************************
# create object with non-numeric prices in lists
error_triggered = False
error_raised = False
try:
_ = ResourcePrice(prices=[7, "4", 2], volumes=[3, 4, 5])
except TypeError:
error_triggered = True
assert error_triggered
error_raised = True
assert error_raised
# *********************************************************************
# create object with non-numeric volumes in lists
error_triggered = False
error_raised = False
try:
_ = ResourcePrice(prices=[7, 3, 2], volumes=[4, "3", 2])
except TypeError:
error_triggered = True
assert error_triggered
error_raised = True
assert error_raised
# *********************************************************************
# create object with mismatched price and volume lists
error_triggered = False
error_raised = False
try:
_ = ResourcePrice(prices=[7, 3, 2], volumes=[5, 7])
except ValueError:
error_triggered = True
assert error_triggered
error_raised = True
assert error_raised
# *********************************************************************
# create object with a price list as an input and an unsupported type
error_triggered = False
error_raised = False
try:
_ = ResourcePrice(prices=[7, 3, 2], volumes="hello")
except TypeError:
error_triggered = True
assert error_triggered
error_raised = True
assert error_raised
# *********************************************************************
# create object with negative prices in lists (no volumes are provided)
error_triggered = False
error_raised = False
try:
_ = ResourcePrice(prices=[7, 3, -2], volumes=None)
except TypeError:
error_triggered = True
assert error_triggered
error_raised = True
assert error_raised
# *********************************************************************
# create object with non-numeric prices in lists (no volumes are provided)
error_triggered = False
error_raised = False
try:
_ = ResourcePrice(prices=[7, 3, "a"], volumes=None)
except TypeError:
error_triggered = True
assert error_triggered
error_raised = True
assert error_raised
# *********************************************************************
# create object with non-numeric prices in lists (no volumes are provided)
error_triggered = False
error_raised = False
try:
_ = ResourcePrice(prices=5, volumes=[7, 3, 4])
except TypeError:
error_triggered = True
assert error_triggered
error_raised = True
assert error_raised
# *********************************************************************
# create object with negative prices
error_triggered = False
error_raised = False
try:
_ = ResourcePrice(prices=-3, volumes=None)
except ValueError:
error_triggered = True
assert error_triggered
error_raised = True
assert error_raised
# *********************************************************************
......