From f7d386caa709744e4d1249e945b7d8bb998abf9a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Pedro=20L=2E=20Magalh=C3=A3es?=
 <pedro.magalhaes@uni-bremen.de>
Date: Fri, 8 Mar 2024 20:30:38 +0100
Subject: [PATCH] Added tests for convex price functions.

---
 src/topupopt/problems/esipp/problem.py |  32 +-
 tests/test_esipp_problem.py            | 443 ++++++++++++++++++++++++-
 tests/test_esipp_resource.py           |   2 +-
 3 files changed, 450 insertions(+), 27 deletions(-)

diff --git a/src/topupopt/problems/esipp/problem.py b/src/topupopt/problems/esipp/problem.py
index 600195c..92d08ca 100644
--- a/src/topupopt/problems/esipp/problem.py
+++ b/src/topupopt/problems/esipp/problem.py
@@ -1,29 +1,19 @@
 # imports
 
 # standard libraries
-
 from numbers import Real
-
 from statistics import mean
-
 import math
 
 # local libraries, external
-
 import pyomo.environ as pyo
 
 # local libraries, internal
-
 from .model import create_model
-
 from ...solvers.interface import SolverInterface
-
 from ...data.finance.invest import discount_factor
-
 from .network import Network, Arcs
-
 from .system import EnergySystem
-
 from .resource import ResourcePrice
            
 # *****************************************************************************
@@ -2041,7 +2031,6 @@ class InfrastructurePlanningProblem(EnergySystem):
         # *********************************************************************
             
         # set of price segments
-        
         set_S = {
             (g, l, q, p, k): tuple(
                 s 
@@ -2064,8 +2053,11 @@ class InfrastructurePlanningProblem(EnergySystem):
             }
                     
         # set of GLKS tuples
-        
-        set_GLQPKS = tuple((*glqpk,s) for glqpk, s in set_S.items())
+        set_GLQPKS = tuple(
+            (*glqpk,s) 
+            for glqpk, s_tuple in set_S.items()
+            for s in s_tuple
+            )
         
         set_GLQPKS_exp = tuple(glqpks for glqpks in set_GLQPKS
                                if glqpks[1] in set_L_exp[glqpks[0]])
@@ -4253,6 +4245,20 @@ def is_peak_total_problem(problem: InfrastructurePlanningProblem) -> bool:
                     return False
                 # if the entries are time invariant, checking one will do
                 break
+        # check export nodes
+        for exp_node_key in net.export_nodes:
+            # is an import node, check if it is time invariant
+            if not net.nodes[exp_node_key][
+                    Network.KEY_NODE_PRICES_TIME_INVARIANT]:
+                return False # is not time invariant
+            # it is time invariant, but is it volume invariant? check any qpk
+            for qpk in net.nodes[exp_node_key][Network.KEY_NODE_PRICES]:
+                if not net.nodes[exp_node_key][
+                        Network.KEY_NODE_PRICES][qpk].is_volume_invariant():
+                    # it is not volume invariant
+                    return False
+                # if the entries are time invariant, checking one will do
+                break
     
     # # check #4: none of the arcs can have proportional losses
     # for key, net in problem.networks.items():
diff --git a/tests/test_esipp_problem.py b/tests/test_esipp_problem.py
index 862e5cf..1b50747 100644
--- a/tests/test_esipp_problem.py
+++ b/tests/test_esipp_problem.py
@@ -1,31 +1,20 @@
 # imports
 
 # standard
-
 import math
-
 from statistics import mean
 
 # local
-
 # import numpy as np
-
 # import networkx as nx
-
 import pyomo.environ as pyo
     
 # import src.topupopt.problems.esipp.utils as utils
-
 from src.topupopt.data.misc.utils import generate_pseudo_unique_key
-
 from src.topupopt.problems.esipp.problem import InfrastructurePlanningProblem
-
 from src.topupopt.problems.esipp.network import Arcs, Network
-
 from src.topupopt.problems.esipp.resource import ResourcePrice
-
 from src.topupopt.problems.esipp.problem import simplify_peak_total_problem
-
 from src.topupopt.problems.esipp.problem import is_peak_total_problem
 
 # *****************************************************************************
@@ -404,7 +393,7 @@ class TestESIPPProblem:
         
         # validation
         
-        # the arc should be installed since it is the only feasible solution
+        # 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
         
@@ -559,7 +548,7 @@ class TestESIPPProblem:
             
         # validation
         
-        # the arc should be installed since it is the only feasible solution
+        # 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
         
@@ -626,10 +615,438 @@ class TestESIPPProblem:
         #     print(12)
         #     12+3
         # print('Got stdout: "{0}"'.format(f.getvalue()))
+    
+    # *************************************************************************
+    # *************************************************************************
+    
+    def test_problem_increasing_imp_prices(self):
+        
+        # scenario
+        q = 0
+        # time
+        number_intervals = 1
+        # periods
+        number_periods = 1
+    
+        # 2 nodes: one import, one regular
+        mynet = Network()
+        
+        # import node        
+        node_IMP = generate_pseudo_unique_key(mynet.nodes())
+        mynet.add_import_node(
+            node_key=node_IMP, 
+            prices={
+                (q,p,k): ResourcePrice(
+                    prices=[1.0, 2.0],
+                    volumes=[0.5, None]
+                    )
+                for p in range(number_periods)
+                for k in range(number_intervals)
+                }
+            )
+        
+        # other nodes
+        node_A = generate_pseudo_unique_key(mynet.nodes())
+        mynet.add_source_sink_node(
+            node_key=node_A, 
+            base_flow={(q, 0): 1.0} 
+            )
+        
+        # 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)
         
+        # identify node types
+        mynet.identify_node_types()
+        
+        # no sos, regular time intervals
+        ipp = self.build_solve_ipp(
+            # solver=solver,
+            solver_options={},
+            # use_sos_arcs=use_sos_arcs,
+            # arc_sos_weight_key=sos_weight_key,
+            # arc_use_real_variables_if_possible=use_real_variables_if_possible,
+            # use_sos_sense=use_sos_sense,
+            # sense_sos_weight_key=sense_sos_weight_key,
+            # sense_use_real_variables_if_possible=sense_use_real_variables_if_possible,
+            # sense_use_arc_interfaces=use_arc_interfaces,
+            perform_analysis=False,
+            plot_results=False, # True,
+            print_solver_output=False,
+            # irregular_time_intervals=irregular_time_intervals,
+            networks={'mynet': mynet},
+            number_intraperiod_time_intervals=number_intervals,
+            static_losses_mode=True, # just to reach a line,
+            mandatory_arcs=[],
+            max_number_parallel_arcs={},
+            # init_aux_sets=init_aux_sets,
+            simplify_problem=False,
+            reporting_periods={0: (0,)},
+            discount_rates={0: (0.0,)}
+            )
+        
+        assert not is_peak_total_problem(ipp)
+        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
         
         # *********************************************************************
         # *********************************************************************
+        
+        # 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 1.0, 0.0 and 2.0
+        assert math.isclose(
+            pyo.value(
+                ipp.instance.var_v_glljqk[('mynet', node_IMP, node_A, 0, q, 0)]
+                ),
+            2.0,
+            abs_tol=1e-6)
+        
+        # arc amplitude should be two
+        assert math.isclose(
+            pyo.value(
+                ipp.instance.var_v_amp_gllj[('mynet', node_IMP, node_A, 0)]
+                ),
+            2.0,
+            abs_tol=0.01)
+        
+        # capex should be four
+        assert math.isclose(pyo.value(ipp.instance.var_capex), 4.0, abs_tol=1e-3)
+        
+        # sdncf should be -3.5
+        assert math.isclose(
+            pyo.value(ipp.instance.var_sdncf_q[q]), -3.5, abs_tol=1e-3
+            )
+        
+        # the objective function should be -7.5
+        assert math.isclose(pyo.value(ipp.instance.obj_f), -7.5, abs_tol=1e-3)
+    
+    # *************************************************************************
+    # *************************************************************************
+    
+    def test_problem_decreasing_exp_prices(self):
+        
+        # scenario
+        q = 0
+        # time
+        number_intervals = 1
+        # periods
+        number_periods = 1
+    
+        # 2 nodes: one export, one regular
+        mynet = Network()
+        
+        # import node        
+        node_EXP = generate_pseudo_unique_key(mynet.nodes())
+        mynet.add_export_node(
+            node_key=node_EXP, 
+            prices={
+                (q,p,k): ResourcePrice(
+                    prices=[2.0, 1.0],
+                    volumes=[0.5, None]
+                    )
+                for p in range(number_periods)
+                for k in range(number_intervals)
+                }
+            )
+        
+        # other nodes
+        node_A = generate_pseudo_unique_key(mynet.nodes())
+        mynet.add_source_sink_node(
+            node_key=node_A, 
+            base_flow={(q, 0): -1.0} 
+            )
+        
+        # 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_A, 
+            node_key_b=node_EXP,
+            arcs=arc_tech_IA)
+        
+        # identify node types
+        mynet.identify_node_types()
+        
+        # no sos, regular time intervals
+        ipp = self.build_solve_ipp(
+            # solver=solver,
+            solver_options={},
+            # use_sos_arcs=use_sos_arcs,
+            # arc_sos_weight_key=sos_weight_key,
+            # arc_use_real_variables_if_possible=use_real_variables_if_possible,
+            # use_sos_sense=use_sos_sense,
+            # sense_sos_weight_key=sense_sos_weight_key,
+            # sense_use_real_variables_if_possible=sense_use_real_variables_if_possible,
+            # sense_use_arc_interfaces=use_arc_interfaces,
+            perform_analysis=False,
+            plot_results=False, # True,
+            print_solver_output=False,
+            # irregular_time_intervals=irregular_time_intervals,
+            networks={'mynet': mynet},
+            number_intraperiod_time_intervals=number_intervals,
+            static_losses_mode=True, # just to reach a line,
+            mandatory_arcs=[],
+            max_number_parallel_arcs={},
+            # init_aux_sets=init_aux_sets,
+            simplify_problem=False,
+            reporting_periods={0: (0,)},
+            discount_rates={0: (0.0,)}
+            )
+        
+        assert not is_peak_total_problem(ipp)
+        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
+        
+        # *********************************************************************
+        # *********************************************************************
+        
+        # validation
+        
+        # the arc should be installed since it is required for feasibility
+        assert True in ipp.networks['mynet'].edges[(node_A, node_EXP, 0)][
+            Network.KEY_ARC_TECH].options_selected
+        
+        # the flows should be 1.0, 0.0 and 2.0
+        assert math.isclose(
+            pyo.value(
+                ipp.instance.var_v_glljqk[('mynet', node_A, node_EXP, 0, q, 0)]
+                ),
+            1.0,
+            abs_tol=1e-6)
+        
+        # arc amplitude should be two
+        assert math.isclose(
+            pyo.value(
+                ipp.instance.var_v_amp_gllj[('mynet', node_A, node_EXP, 0)]
+                ),
+            1.0,
+            abs_tol=0.01)
+        
+        # capex should be four
+        assert math.isclose(pyo.value(ipp.instance.var_capex), 3.0, abs_tol=1e-3)
+        
+        # sdncf should be 1.0
+        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), -2.0, abs_tol=1e-3)
+    
+    # *************************************************************************
+    # *************************************************************************
+    
+    def test_problem_increasing_imp_decreasing_exp_prices(self):
+        
+        # scenario
+        q = 0
+        # time
+        number_intervals = 2
+        # periods
+        number_periods = 1
+    
+        # 3 nodes: one import, one export, one regular
+        mynet = Network()
+        
+        # import node        
+        node_IMP = generate_pseudo_unique_key(mynet.nodes())
+        mynet.add_import_node(
+            node_key=node_IMP, 
+            prices={
+                (q,p,k): ResourcePrice(
+                    prices=[1.0, 2.0],
+                    volumes=[0.5, None]
+                    )
+                for p in range(number_periods)
+                for k in range(number_intervals)
+                }
+            )
+        
+        # export node
+        node_EXP = generate_pseudo_unique_key(mynet.nodes())
+        mynet.add_export_node(
+            node_key=node_EXP, 
+            prices={
+                (q,p,k): ResourcePrice(
+                    prices=[2.0, 1.0],
+                    volumes=[0.5, None]
+                    )
+                for p in range(number_periods)
+                for k in range(number_intervals)
+                }
+            )
+        
+        # other nodes
+        node_A = generate_pseudo_unique_key(mynet.nodes())
+        mynet.add_source_sink_node(
+            node_key=node_A, 
+            base_flow={(q, 0): 1.0, (q, 1): -1.0} 
+            )
+        
+        # arc IA
+        arc_tech_IA = Arcs(
+            name='any',
+            efficiency={(q, 0): 0.5, (q, 1): 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)
+        
+        # arc AE
+        arc_tech_AE = Arcs(
+            name='any',
+            efficiency={(q, 0): 0.5, (q, 1): 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_A, 
+            node_key_b=node_EXP,
+            arcs=arc_tech_AE)
+        
+        # identify node types
+        mynet.identify_node_types()
+        
+        # no sos, regular time intervals
+        ipp = self.build_solve_ipp(
+            # solver=solver,
+            solver_options={},
+            # use_sos_arcs=use_sos_arcs,
+            # arc_sos_weight_key=sos_weight_key,
+            # arc_use_real_variables_if_possible=use_real_variables_if_possible,
+            # use_sos_sense=use_sos_sense,
+            # sense_sos_weight_key=sense_sos_weight_key,
+            # sense_use_real_variables_if_possible=sense_use_real_variables_if_possible,
+            # sense_use_arc_interfaces=use_arc_interfaces,
+            perform_analysis=False,
+            plot_results=False, # True,
+            print_solver_output=False,
+            # irregular_time_intervals=irregular_time_intervals,
+            networks={'mynet': mynet},
+            number_intraperiod_time_intervals=number_intervals,
+            static_losses_mode=True, # just to reach a line,
+            mandatory_arcs=[],
+            max_number_parallel_arcs={},
+            # init_aux_sets=init_aux_sets,
+            simplify_problem=False,
+            reporting_periods={0: (0,)},
+            discount_rates={0: (0.0,)}
+            )
+        
+        assert not is_peak_total_problem(ipp)
+        assert ipp.results['Problem'][0]['Number of constraints'] == 23
+        assert ipp.results['Problem'][0]['Number of variables'] == 26
+        assert ipp.results['Problem'][0]['Number of nonzeros'] == 57
+        
+        # *********************************************************************
+        # *********************************************************************
+        
+        # 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 arc should be installed since it is required for feasibility
+        assert True in ipp.networks['mynet'].edges[(node_A, node_EXP, 0)][
+            Network.KEY_ARC_TECH].options_selected
+        
+        # interval 0: import only
+        assert math.isclose(
+            pyo.value(
+                ipp.instance.var_v_glljqk[('mynet', node_IMP, node_A, 0, q, 0)]
+                ),
+            2.0,
+            abs_tol=1e-6
+            )
+        assert math.isclose(
+            pyo.value(
+                ipp.instance.var_v_glljqk[('mynet', node_A, node_EXP, 0, q, 0)]
+                ),
+            0.0,
+            abs_tol=1e-6
+            )
+        # interval 1: export only
+        assert math.isclose(
+            pyo.value(
+                ipp.instance.var_v_glljqk[('mynet', node_IMP, node_A, 0, q, 1)]
+                ),
+            0.0,
+            abs_tol=1e-6
+            )
+        assert math.isclose(
+            pyo.value(
+                ipp.instance.var_v_glljqk[('mynet', node_A, node_EXP, 0, q, 1)]
+                ),
+            1.0,
+            abs_tol=1e-6
+            )
+        
+        # IA amplitude
+        assert math.isclose(
+            pyo.value(
+                ipp.instance.var_v_amp_gllj[('mynet', node_IMP, node_A, 0)]
+                ),
+            2.0,
+            abs_tol=0.01)
+        # AE amplitude
+        assert math.isclose(
+            pyo.value(
+                ipp.instance.var_v_amp_gllj[('mynet', node_A, node_EXP, 0)]
+                ),
+            1.0,
+            abs_tol=0.01)
+        
+        # capex should be 7.0: 4+3
+        assert math.isclose(pyo.value(ipp.instance.var_capex), 7.0, abs_tol=1e-3)
+        
+        # sdncf should be -2.5: -3.5+1.0
+        assert math.isclose(
+            pyo.value(ipp.instance.var_sdncf_q[q]), -2.5, abs_tol=1e-3
+            )
+        
+        # the objective function should be -9.5: -7.5-2.5
+        assert math.isclose(pyo.value(ipp.instance.obj_f), -9.5, abs_tol=1e-3)
+    
+    # *************************************************************************
+    # *************************************************************************
     
 # *****************************************************************************
 # *****************************************************************************
\ No newline at end of file
diff --git a/tests/test_esipp_resource.py b/tests/test_esipp_resource.py
index 642a3e1..8685e8b 100644
--- a/tests/test_esipp_resource.py
+++ b/tests/test_esipp_resource.py
@@ -7,7 +7,7 @@ from src.topupopt.problems.esipp.resource import are_prices_time_invariant
 # *****************************************************************************
 
 class TestResourcePrice:
-    
+        
     # *************************************************************************
     # *************************************************************************
     
-- 
GitLab