diff --git a/src/topupopt/problems/esipp/blocks/delta.py b/src/topupopt/problems/esipp/blocks/delta.py index 59b418feb43902d31282c8666f7417ce3071b7f1..2cdc79bff68a727a725a8c8db28b3caacdf246d2 100644 --- a/src/topupopt/problems/esipp/blocks/delta.py +++ b/src/topupopt/problems/esipp/blocks/delta.py @@ -25,11 +25,6 @@ def price_delta_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) # ********************************************************************* # ********************************************************************* @@ -66,28 +61,21 @@ def price_delta_block( bounds=bounds_var_trans_flows ) - # import and export flows - def rule_var_segment_usage_s(b, s): - if b.param_price_function_is_convex: - return pyo.UnitInterval - else: - return pyo.Binary + # segment usage variables b.var_segment_usage_s = pyo.Var( b.set_S, - within=rule_var_segment_usage_s, + within=pyo.UnitInterval, ) # ********************************************************************* # ********************************************************************* - # import flow costs and export flow revenues + # 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] - if s == 0 else - b.var_segment_usage_s[s]*(b.param_p_s[s]*b.param_v_max_s[s]-b.param_p_s[s-1]*b.param_v_max_s[s-1])) + 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)] @@ -95,9 +83,7 @@ def price_delta_block( else: return ( sum( - (b.var_segment_usage_s[s]*b.param_p_s[s]*b.param_v_max_s[s] - if s == 0 else - b.var_segment_usage_s[s]*(b.param_p_s[s]*b.param_v_max_s[s]-b.param_p_s[s-1]*b.param_v_max_s[s-1])) + 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)] @@ -106,7 +92,7 @@ def price_delta_block( rule=rule_constr_trans_monetary_flows ) - # imported and exported flows: volumes + # 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] @@ -114,7 +100,7 @@ def price_delta_block( ) == b.var_trans_flows b.constr_trans_flows_seg = pyo.Constraint(rule=rule_constr_trans_flows_seg) - # imported and exported flows: energy system + # 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( @@ -140,22 +126,64 @@ def price_delta_block( if not b.param_price_function_is_convex: # this means the segment usage variables are binary - raise NotImplementedError + # 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 + ) - # # 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: - # # last segment, skip - # # convex, skip - # return pyo.Constraint.Skip - # return ( - # b.var_segment_usage_s[s] >= - # b.var_segment_usage_s[s+1] - # ) - # b.constr_delta_summing_logic = pyo.Constraint( - # b.set_S, rule=rule_constr_delta_summing_logic - # ) # ********************************************************************* # ********************************************************************* @@ -186,8 +214,6 @@ def price_delta_no_block( enable_initialisation: bool = True ): - raise NotImplementedError - # auxiliary set for pyomo model.set_GLQPK = model.set_GL_exp_imp*model.set_QPK @@ -241,128 +267,146 @@ def price_delta_no_block( # ************************************************************************* # import and export flows - def bounds_var_trans_flows_glqpks(m, g, l, q, p, k, s): - # return (0, m.param_v_max_glqpks[(g, l, q, p, k, s)]) - if (g, l, q, p, k, s) in m.param_v_max_glqpks: - # predefined finite capacity - return (0, m.param_v_max_glqpks[(g, l, q, p, k, s)]) - else: - # infinite capacity + 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_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, 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 + + + # 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_trans_flows_glqpks[(g, l, q, p, k, s)] - * m.param_p_glqpks[(g, l, q, p, k, s)] + 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)] + == 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)] + 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)] + == 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): + 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)] + 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_trans_flows_sys = pyo.Constraint( + model.set_GLQPK, + rule=rule_constr_trans_flows_sys + ) # ************************************************************************* # ************************************************************************* # 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 + + # variables to indicate if the segments are active + model.var_segment_is_full = 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 + # 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_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)] + m.var_segment_usage_glqpks[(g,l,q,p,k,s)] >= + m.var_segment_is_full[(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 + model.constr_nonconvex_p2 = pyo.Constraint( + model.set_GLQPKS, rule=rule_constr_nonconvex_p2 ) - # 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)]: + 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 - # 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)] + 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_delta_summing_logic = pyo.Constraint( - model.set_GLQPKS, rule=rule_constr_delta_summing_logic + model.constr_nonconvex_p3 = pyo.Constraint( + model.set_GLQPKS, rule=rule_constr_nonconvex_p3 ) - # 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)]: + 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 - # 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)] + m.var_segment_usage_glqpks[(g,l,q,p,k,s+1)] <= + m.var_segment_is_full[(g,l,q,p,k,s+1)] ) - # 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 + model.constr_nonconvex_p4 = pyo.Constraint( + model.set_GLQPKS, rule=rule_constr_nonconvex_p4 ) # ***************************************************************************** diff --git a/src/topupopt/problems/esipp/blocks/prices.py b/src/topupopt/problems/esipp/blocks/prices.py index 4c38b3fe092c8b0d23b56cb0030acc86e52ddd87..f4f946902c9cc7c569fa566a6e19c40a0b902a9f 100644 --- a/src/topupopt/problems/esipp/blocks/prices.py +++ b/src/topupopt/problems/esipp/blocks/prices.py @@ -39,7 +39,7 @@ def add_price_functions( if node_price_model == NODE_PRICE_LAMBDA: raise NotImplementedError elif node_price_model == NODE_PRICE_DELTA: - raise NotImplementedError + price_delta_no_block(model, **kwargs) else: price_other_no_block(model, **kwargs) diff --git a/src/topupopt/problems/esipp/utils.py b/src/topupopt/problems/esipp/utils.py index a61b35a9ac94ad5ae7228308ab8baeb8703d7542..bffc0b94be97a0c9892c9a3420fb92e2a8eece82 100644 --- a/src/topupopt/problems/esipp/utils.py +++ b/src/topupopt/problems/esipp/utils.py @@ -152,35 +152,60 @@ def statistics(ipp: InfrastructurePlanningProblem, } else: # not in a block - # 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)] + if ipp.node_price_model == NODE_PRICE_DELTA: + # 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] ) - *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)] + 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] ) - *ipp.instance.param_c_time_qpk[qpk] - ) - for qpk in ipp.time_frame.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] diff --git a/tests/test_esipp_prices.py b/tests/test_esipp_prices.py index 8cd73ddcfb2e7279fb250e8b51975c92ba3dfd26..2652bd64451e91152af4a5861a561a0042e524b8 100644 --- a/tests/test_esipp_prices.py +++ b/tests/test_esipp_prices.py @@ -25,6 +25,9 @@ from src.topupopt.problems.esipp.blocks.prices import NODE_PRICE_OTHER, NODE_PRI # ***************************************************************************** # ***************************************************************************** +# TODO: test time-varying tariffs (convex/non-convex) +# TODO: check problem sizes + class TestESIPPProblem: # ************************************************************************* @@ -102,14 +105,25 @@ class TestESIPPProblem: ) assert not ipp.has_peak_total_assessments() - - print('hey') - print(ipp.results["Problem"][0]) - # if node_price_model != NODE_PRICE_OTHER: - # ipp.instance.pprint() - # assert ipp.results["Problem"][0]["Number of constraints"] == 24 - # assert ipp.results["Problem"][0]["Number of variables"] == 22 - # assert ipp.results["Problem"][0]["Number of nonzeros"] == 49 + # print('hey') + # print((use_prices_block, node_price_model)) + # print(ipp.results["Problem"][0]) + if (use_prices_block, node_price_model) == (True, NODE_PRICE_OTHER): + assert ipp.results["Problem"][0]["Number of constraints"] == 10 + assert ipp.results["Problem"][0]["Number of variables"] == 11 + assert ipp.results["Problem"][0]["Number of nonzeros"] == 20 + elif (use_prices_block, node_price_model) == (False, NODE_PRICE_OTHER): + assert ipp.results["Problem"][0]["Number of constraints"] == 10 + assert ipp.results["Problem"][0]["Number of variables"] == 11 + assert ipp.results["Problem"][0]["Number of nonzeros"] == 20 + elif (use_prices_block, node_price_model) == (True, NODE_PRICE_DELTA): + assert ipp.results["Problem"][0]["Number of constraints"] == 11 + assert ipp.results["Problem"][0]["Number of variables"] == 12 + assert ipp.results["Problem"][0]["Number of nonzeros"] == 22 + elif (use_prices_block, node_price_model) == (False, NODE_PRICE_DELTA): + assert ipp.results["Problem"][0]["Number of constraints"] == 15 + assert ipp.results["Problem"][0]["Number of variables"] == 14 + assert ipp.results["Problem"][0]["Number of nonzeros"] == 30 # ********************************************************************* # ********************************************************************* @@ -149,7 +163,8 @@ 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), @@ -157,7 +172,8 @@ class TestESIPPProblem: (True, NODE_PRICE_LAMBDA), (False, NODE_PRICE_OTHER), (False, NODE_PRICE_DELTA), - (False, NODE_PRICE_LAMBDA)] + (False, NODE_PRICE_LAMBDA) + ] ) def test_problem_decreasing_imp_prices(self, use_prices_block, node_price_model): @@ -220,6 +236,25 @@ class TestESIPPProblem: ) assert not ipp.has_peak_total_assessments() + # print('hey') + # print((use_prices_block, node_price_model)) + # print(ipp.results["Problem"][0]) + if (use_prices_block, node_price_model) == (True, NODE_PRICE_OTHER): + assert ipp.results["Problem"][0]["Number of constraints"] == 14 + assert ipp.results["Problem"][0]["Number of variables"] == 13 + assert ipp.results["Problem"][0]["Number of nonzeros"] == 28 + elif (use_prices_block, node_price_model) == (False, NODE_PRICE_OTHER): + assert ipp.results["Problem"][0]["Number of constraints"] == 14 + assert ipp.results["Problem"][0]["Number of variables"] == 13 + assert ipp.results["Problem"][0]["Number of nonzeros"] == 28 + elif (use_prices_block, node_price_model) == (True, NODE_PRICE_DELTA): + assert ipp.results["Problem"][0]["Number of constraints"] == 15 + assert ipp.results["Problem"][0]["Number of variables"] == 14 + assert ipp.results["Problem"][0]["Number of nonzeros"] == 30 + elif (use_prices_block, node_price_model) == (False, NODE_PRICE_DELTA): + assert ipp.results["Problem"][0]["Number of constraints"] == 15 + assert ipp.results["Problem"][0]["Number of variables"] == 14 + assert ipp.results["Problem"][0]["Number of nonzeros"] == 30 # ********************************************************************* # ********************************************************************* @@ -256,6 +291,140 @@ class TestESIPPProblem: # the objective function should be -7.5 assert math.isclose(pyo.value(ipp.instance.obj_f), -6.5, abs_tol=1e-3) + + # ************************************************************************* + # ************************************************************************* + + @pytest.mark.parametrize( + "use_prices_block, node_price_model", + [(True, NODE_PRICE_OTHER), + (True, NODE_PRICE_DELTA), + (True, NODE_PRICE_LAMBDA), + (False, NODE_PRICE_OTHER), + (False, NODE_PRICE_DELTA), + (False, NODE_PRICE_LAMBDA) + ] + ) + def test_problem_decreasing_imp_prices2(self, use_prices_block, node_price_model): + + # assessment + q = 0 + + tf = EconomicTimeFrame( + discount_rate=0.0, + reporting_periods={q: (0,)}, + reporting_period_durations={q: (365 * 24 * 3600,)}, + time_intervals={q: (0,)}, + time_interval_durations={q: (1,)}, + ) + + # 2 nodes: one import, one regular + mynet = Network() + + # import node + node_IMP = 'I' + mynet.add_import_node( + node_key=node_IMP, + prices={ + qpk: ResourcePrice(prices=[2.0, 1.0], volumes=[0.5, 3.0]) + for qpk in tf.qpk() + }, + ) + + # other nodes + node_A = 'A' + mynet.add_source_sink_node(node_key=node_A, base_flow={(q, 0): 0.25}) + + # arc IA + arc_tech_IA = Arcs( + name="any", + efficiency={(q, 0): 0.5}, + efficiency_reverse=None, + static_loss=None, + capacity=[3], + minimum_cost=[2], + specific_capacity_cost=1, + capacity_is_instantaneous=False, + validate=False, + ) + mynet.add_directed_arc(node_key_a=node_IMP, node_key_b=node_A, arcs=arc_tech_IA) + + # no sos, regular time intervals + ipp = build_solve_ipp( + solver_options={}, + perform_analysis=False, + plot_results=False, # True, + print_solver_output=False, + time_frame=tf, + networks={"mynet": mynet}, + static_losses_mode=True, # just to reach a line, + mandatory_arcs=[], + max_number_parallel_arcs={}, + simplify_problem=False, + use_prices_block=use_prices_block, + node_price_model=node_price_model + ) + + assert not ipp.has_peak_total_assessments() + # print('hey') + # # print((use_prices_block, node_price_model)) + # print(ipp.results["Problem"][0]) + # # capex should be four + # print(pyo.value(ipp.instance.var_capex)) + # print(pyo.value(ipp.instance.var_sdncf_q[q])) + # print(pyo.value(ipp.instance.obj_f)) + if (use_prices_block, node_price_model) == (True, NODE_PRICE_OTHER): + assert ipp.results["Problem"][0]["Number of constraints"] == 14 + assert ipp.results["Problem"][0]["Number of variables"] == 13 + assert ipp.results["Problem"][0]["Number of nonzeros"] == 28 + elif (use_prices_block, node_price_model) == (False, NODE_PRICE_OTHER): + assert ipp.results["Problem"][0]["Number of constraints"] == 14 + assert ipp.results["Problem"][0]["Number of variables"] == 13 + assert ipp.results["Problem"][0]["Number of nonzeros"] == 28 + elif (use_prices_block, node_price_model) == (True, NODE_PRICE_DELTA): + assert ipp.results["Problem"][0]["Number of constraints"] == 15 + assert ipp.results["Problem"][0]["Number of variables"] == 14 + assert ipp.results["Problem"][0]["Number of nonzeros"] == 30 + elif (use_prices_block, node_price_model) == (False, NODE_PRICE_DELTA): + assert ipp.results["Problem"][0]["Number of constraints"] == 15 + assert ipp.results["Problem"][0]["Number of variables"] == 14 + assert ipp.results["Problem"][0]["Number of nonzeros"] == 30 + + # ********************************************************************* + # ********************************************************************* + + # validation + + # the arc should be installed since it is required for feasibility + assert ( + True + in ipp.networks["mynet"] + .edges[(node_IMP, node_A, 0)][Network.KEY_ARC_TECH] + .options_selected + ) + + # the flows should be 0.5 + assert math.isclose( + pyo.value(ipp.instance.var_v_glljqk[("mynet", node_IMP, node_A, 0, q, 0)]), + 0.5, + abs_tol=1e-6, + ) + + # arc amplitude should be 0.5 + assert math.isclose( + pyo.value(ipp.instance.var_v_amp_gllj[("mynet", node_IMP, node_A, 0)]), + 0.5, + abs_tol=0.01, + ) + + # capex should be four + assert math.isclose(pyo.value(ipp.instance.var_capex), 2.5, abs_tol=1e-3) + + # sdncf should be -2.5 + assert math.isclose(pyo.value(ipp.instance.var_sdncf_q[q]), -1.0, abs_tol=1e-3) + + # the objective function should be -7.5 + assert math.isclose(pyo.value(ipp.instance.obj_f), -3.5, abs_tol=1e-3) # ************************************************************************* # *************************************************************************