From fe4f7bc0a0587f8798234e3319eb5ffba479c3b7 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 02:17:58 +0100
Subject: [PATCH] Revised classes for investments in generic arcs and in
district heating trenches. Updated converter and dynsys tests.
---
src/topupopt/data/dhn/network.py | 315 ++-
src/topupopt/data/dhn/utils.py | 68 +-
src/topupopt/data/finance/invest.py | 3 +
src/topupopt/problems/esipp/network.py | 33 +-
src/topupopt/problems/esipp/problem.py | 6 +-
tests/examples_dhn.py | 1564 -----------
tests/examples_dynsys.py | 2366 -----------------
tests/test_all.py | 30 -
tests/test_converter.py | 292 --
tests/test_data_finance.py | 207 +-
tests/test_dhn.py | 168 +-
tests/test_dhn_utils.py | 227 +-
...s_converter.py => test_esipp_converter.py} | 48 +-
tests/test_esipp_dynsys.py | 2231 ++++++++++++++++
tests/test_esipp_network.py | 2061 ++++++++++++++
tests/test_esipp_problem.py | 30 +-
16 files changed, 5215 insertions(+), 4434 deletions(-)
delete mode 100644 tests/examples_dhn.py
delete mode 100644 tests/examples_dynsys.py
delete mode 100644 tests/test_converter.py
rename tests/{examples_converter.py => test_esipp_converter.py} (88%)
create mode 100644 tests/test_esipp_dynsys.py
create mode 100644 tests/test_esipp_network.py
diff --git a/src/topupopt/data/dhn/network.py b/src/topupopt/data/dhn/network.py
index 3be56ef..a8839ce 100644
--- a/src/topupopt/data/dhn/network.py
+++ b/src/topupopt/data/dhn/network.py
@@ -10,9 +10,12 @@ from numbers import Real
from topupheat.pipes.trenches import SupplyReturnPipeTrench
# local, internal imports
-from ...data.finance.invest import Investment
from ...problems.esipp.network import ArcsWithoutProportionalLosses
+# TODO: add way to define polarity of cash flows
+
+from ...data.finance.utils import ArcInvestments
+
#******************************************************************************
#******************************************************************************
@@ -31,6 +34,7 @@ KEY_HHT_STD_PIPES = 'pipe_tuple'
# *****************************************************************************
class PipeTrenchOptions(ArcsWithoutProportionalLosses):
+ "A class for defining investments in district heating trenches."
def __init__(
self,
@@ -203,150 +207,201 @@ class PipeTrenchOptions(ArcsWithoutProportionalLosses):
self.static_loss = {
(0, scenario_key, 0): hts
}
-
-# *****************************************************************************
-# *****************************************************************************
-
-class ExistingPipeTrench(PipeTrenchOptions):
-
- def __init__(self, option_selected: int, **kwargs):
-
- PipeTrenchOptions.__init__(
- self,
- minimum_cost=[0 for i in range(kwargs['trench'].number_options())],
- **kwargs)
- # define the option that already exists
- self.options_selected[option_selected] = True
# *****************************************************************************
# *****************************************************************************
-# TODO: add way to define polarity of cash flows
-
-class PipeTrenchInvestmentOptions(PipeTrenchOptions):
+class PipeTrenchInvestments(ArcInvestments, PipeTrenchOptions):
+ "A class for defining investments in district heating trenches."
def __init__(
- self,
- discount_rates: list, # 1 per period
- # extra
- investment_besides_pipes: float or list, #
- investment_period: float or list,
- investment_longevity: float or list,
- commissioning_delay_after_investment: float or list,
- fixed_operational_cash_flows: float or list = None,
- start_operational_cash_flows: float or list = None,
- salvage_value_method: str = 'other',
- **kwargs):
- # prepare the instance as usual: minimum cost = cost of pipes
- PipeTrenchOptions.__init__(
self,
+ trench: SupplyReturnPipeTrench,
+ name: str,
+ length: float,
+ investments: tuple,
+ static_loss: dict = None,
+ specific_capacity_cost: float or list = None,
+ capacity_is_instantaneous: bool = False,
+ unit_conversion_factor: float = 1.0,
**kwargs
+ ):
+
+ # store the unit conversion
+ self.unit_conversion_factor = unit_conversion_factor
+ # keep the trench object
+ self.trench = trench
+ # keep the trench length
+ self.length = (
+ [length for i in range(trench.number_options())]
+ if trench.vector_mode else
+ length
)
- # investments
- self.invs = [
- Investment(discount_rates=discount_rates)
- for i in range(self.number_options())
- ]
- # add cost of pipes (self.minimum cost) and other items
- self.add_investment(
- investment=list(
- map(sum,zip(self.minimum_cost, investment_besides_pipes))
+ # determine the rated heat capacity
+ rhc = trench.rated_heat_capacity(
+ unit_conversion_factor=unit_conversion_factor
+ )
+ # initialise the object using the mother class
+ ArcInvestments.__init__(
+ self,
+ investments=investments,
+ name=name,
+ efficiency=None,
+ efficiency_reverse=None,
+ static_loss=static_loss,
+ capacity=[rhc] if isinstance(rhc, Real) else rhc,
+ specific_capacity_cost=(
+ 0
+ if type(specific_capacity_cost) == type(None) else
+ specific_capacity_cost
),
- investment_period=investment_period,
- investment_longevity=investment_longevity,
- commissioning_delay_after_investment=(
- commissioning_delay_after_investment),
- salvage_value_method=salvage_value_method
+ capacity_is_instantaneous=False,
+ validate=False
)
- # define the minimum cost as the net present value of the investments
- self.update_minimum_cost()
-
- # *************************************************************************
- # *************************************************************************
+
+ # # *************************************************************************
+ # # *************************************************************************
- def update_minimum_cost(self):
- """Sets the minimum cost as each investment\'s net prevent value."""
- # define the minimum cost
- self.minimum_cost = tuple(inv.net_present_value() for inv in self.invs)
+ # def set_minimum_cost(self, minimum_cost = None):
- # *************************************************************************
- # *************************************************************************
-
- def add_investment(
- self,
- investment: float or list, #
- investment_period: float or list,
- investment_longevity: float or list,
- commissioning_delay_after_investment: float or list,
- salvage_value_method: str = 'other',
- update_minimum_cost: bool = False
- ):
+ # # minimum arc cost
+ # # if no external minimum cost list was provided, calculate it
+ # if type(minimum_cost) == type(None):
+ # # use the specific pipe cost that features in the database
+ # if self.trench.vector_mode:
+ # # multiple options
+ # self.minimum_cost = tuple(
+ # (pipe.sp*length # twin pipes: one twin pipe
+ # if self.trench.twin_pipes else
+ # pipe.sp*length*2) # single pipes: two single pipes
+ # for pipe, length in zip(
+ # self.trench.supply_pipe,
+ # self.length
+ # )
+ # )
+ # else: # only one option
+ # self.minimum_cost = (self.trench.supply_pipe.sp*self.length,)
+ # else: # use an external minimum cost
+ # self.minimum_cost = tuple(minimum_cost)
+
+ # # *************************************************************************
+ # # *************************************************************************
- # check the inputs
- if (isinstance(investment, Real) and
- isinstance(investment_period, Real) and
- isinstance(investment_longevity, Real) and
- isinstance(commissioning_delay_after_investment, Real)):
- # real numbers as inputs: normal mode
- vector_mode = False
- # store the details
- # self.investment = investment
- # self.investment_period = investment_period
- # self.investment_longevity = investment_longevity
- # self.commissioning_delay_after_investment = (
- # commissioning_delay_after_investment)
- elif (type(investment) == list and
- type(investment_period) == list and
- type(investment_longevity) == list and
- type(commissioning_delay_after_investment) == list):
- # lists as inputs: vector mode
- vector_mode = True
- # # store the details
- # self.investment = investment
- # self.investment_period = investment_period
- # self.investment_longevity = investment_longevity
- # self.commissioning_delay_after_investment = (
- # commissioning_delay_after_investment)
- else:
- raise TypeError('Incompatible inputs.')
- # make sure the investment and trench details are consistent
- if vector_mode != self.trench.vector_mode:
- raise TypeError('Inconsistent input types.')
- if vector_mode and self.number_options() != len(investment):
- raise ValueError('Inconsistent inputs.')
- # add the investment
- if vector_mode:
- # multiple options
- for inv, ii, ip, il, cdai in zip(
- self.invs,
- investment,
- investment_period,
- investment_longevity,
- commissioning_delay_after_investment
- ):
- inv.add_investment(
- investment=ii,
- investment_period=ip,
- investment_longevity=il,
- commissioning_delay_after_investment=cdai,
- salvage_value_method=salvage_value_method
- )
- else:
- # single option
- self.invs[0].add_investment(
- investment=investment,
- investment_period=investment_period,
- investment_longevity=investment_longevity,
- commissioning_delay_after_investment=(
- commissioning_delay_after_investment),
- salvage_value_method=salvage_value_method
- )
- # update the minimum cost if desired
- if update_minimum_cost:
- self.update_minimum_cost()
+ # def set_capacity(self, **kwargs):
+ # # retrieve the rated heat capacity
+ # rhc = self.trench.rated_heat_capacity(**kwargs)
+ # if self.trench.vector_mode:
+ # # multiple options, rhc is already a list of values
+ # self.capacity = rhc
+ # else:
+ # # one option, rhc is one value
+ # self.capacity = (rhc,)
- # *************************************************************************
- # *************************************************************************
+ # # *************************************************************************
+ # # *************************************************************************
+
+ # def set_static_losses(
+ # self,
+ # scenario_key,
+ # ground_thermal_conductivity: float or list,
+ # ground_air_heat_transfer_coefficient: float or list,
+ # time_interval_duration: float or list,
+ # temperature_surroundings: float or list,
+ # length: float or list = None,
+ # unit_conversion_factor: float = None,
+ # **kwargs):
+
+ # hts = self.trench.heat_transfer_surroundings(
+ # ground_thermal_conductivity=ground_thermal_conductivity,
+ # ground_air_heat_transfer_coefficient=(
+ # ground_air_heat_transfer_coefficient),
+ # time_interval_duration=time_interval_duration,
+ # temperature_surroundings=temperature_surroundings,
+ # length=(
+ # self.length
+ # if type(length) == type(None) else
+ # length
+ # ),
+ # unit_conversion_factor=(
+ # self.unit_conversion_factor
+ # if type(unit_conversion_factor) == type(None) else
+ # unit_conversion_factor
+ # ),
+ # **kwargs)
+
+ # if self.trench.vector_mode:
+ # # multiple options: hts is a vector
+ # if (hasattr(self, "static_loss") and
+ # type(self.static_loss) != type(None)):
+ # # update the static loss dictionary
+ # if type(hts[0]) == list:
+ # # multiple time intervals
+ # self.static_loss.update({
+ # (h, scenario_key, k): hts[h][k]
+ # for h, hts_h in enumerate(hts)
+ # for k, hts_hk in enumerate(hts_h)
+ # })
+ # else: # not a list: one time interval
+ # self.static_loss.update({
+ # (h, scenario_key, 0): hts[h]
+ # for h, hts_h in enumerate(hts)
+ # })
+ # else:
+ # # no static loss dictionary, create it
+ # if type(hts[0]) == list:
+ # # multiple time intervals
+ # self.static_loss = {
+ # (h, scenario_key, k): hts[h][k]
+ # for h, hts_h in enumerate(hts)
+ # for k, hts_hk in enumerate(hts_h)
+ # }
+ # else: # not a list: one time interval
+ # self.static_loss = {
+ # (h, scenario_key, 0): hts[h]
+ # for h, hts_h in enumerate(hts)
+ # }
+ # else:
+ # # one option: hts might be a number
+ # if (hasattr(self, "static_loss") and
+ # type(self.static_loss) != type(None)):
+ # # update the static loss dictionary
+ # if not isinstance(hts, Real):
+ # # multiple time intervals
+ # self.static_loss.update({
+ # (0, scenario_key, k): hts[k]
+ # for k, hts_k in enumerate(hts)
+ # })
+ # else: # not a list: one time interval
+ # self.static_loss.update({
+ # (0, scenario_key, 0): hts
+ # })
+ # else:
+ # # no static loss dictionary, create it
+ # if not isinstance(hts, Real):
+ # # multiple time intervals
+ # self.static_loss = {
+ # (0, scenario_key, k): hts_k
+ # for k, hts_k in enumerate(hts)
+ # }
+ # else: # not a list: one time interval
+ # self.static_loss = {
+ # (0, scenario_key, 0): hts
+ # }
+# *****************************************************************************
+# *****************************************************************************
+
+class ExistingPipeTrench(PipeTrenchOptions):
+ "A class for existing pipe trenches."
+
+ def __init__(self, option_selected: int, **kwargs):
+ # initialise
+ PipeTrenchOptions.__init__(
+ self,
+ minimum_cost=[0 for i in range(kwargs['trench'].number_options())],
+ **kwargs)
+ # define the option that already exists
+ self.options_selected[option_selected] = True
+
# *****************************************************************************
# *****************************************************************************
\ No newline at end of file
diff --git a/src/topupopt/data/dhn/utils.py b/src/topupopt/data/dhn/utils.py
index df6aee9..8ae75f3 100644
--- a/src/topupopt/data/dhn/utils.py
+++ b/src/topupopt/data/dhn/utils.py
@@ -13,14 +13,78 @@ import numpy as np
from ...problems.esipp.network import Network
from .network import PipeTrenchOptions
+from topupheat.pipes.trenches import SupplyReturnPipeTrench
+from numbers import Real
+
+# *****************************************************************************
+# *****************************************************************************
+
+def cost_pipes(trench: SupplyReturnPipeTrench,
+ length: float or tuple) -> tuple:
+ """
+ Returns the costs of each trench option for a given trench length.
+
+ Parameters
+ ----------
+ trench : SupplyReturnPipeTrench
+ The object describing the trench options.
+ length : float or tuple
+ The trench length in meters.
+
+ Raises
+ ------
+ ValueError
+ This error is raised if unrecognised inputs are inserted.
+
+ Returns
+ -------
+ tuple
+ A tuple of the pipe costs for each option.
+
+ """
+ # use the specific pipe cost that features in the database
+ if trench.vector_mode:
+ # multiple options
+ if (type(length) == tuple and
+ len(length) == trench.number_options()):
+ # multiple trench lengths
+ return tuple(
+ (pipe.sp*length # twin pipes: one twin pipe
+ if trench.twin_pipes else
+ pipe.sp*length*2) # single pipes: two single pipes
+ for pipe, length in zip(trench.supply_pipe, length)
+ )
+ elif isinstance(length, Real):
+ # one trench length
+ return tuple(
+ (pipe.sp*length # twin pipes: one twin pipe
+ if trench.twin_pipes else
+ pipe.sp*length*2) # single pipes: two single pipes
+ for pipe in trench.supply_pipe
+ )
+ else:
+ raise ValueError('Unrecognised input combination.')
+ elif (not trench.vector_mode and isinstance(length, Real)):
+ # only one option
+ return (trench.supply_pipe.sp*length,)
+ else: # only one option
+ raise ValueError('Unrecognised input combination.')
+
+ # # keep the trench length
+ # self.length = (
+ # [length for i in range(trench.number_options())]
+ # if trench.vector_mode else
+ # length
+ # )
# *****************************************************************************
# *****************************************************************************
-def summarise_network_by_arc_technology(
+def summarise_network_by_pipe_technology(
network: Network,
print_output: bool = False
) -> dict:
+ "A method to summarise a network by pipe technology."
# *************************************************************************
# *************************************************************************
@@ -189,7 +253,7 @@ def plot_network_layout(network: Network,
# add base map
if include_basemap:
- # TODO: reach this statement in tests
+
cx.add_basemap(ax,
zoom=basemap_zoom_level,
source=cx.providers.OpenStreetMap.Mapnik,
diff --git a/src/topupopt/data/finance/invest.py b/src/topupopt/data/finance/invest.py
index ba96df7..ca682b2 100644
--- a/src/topupopt/data/finance/invest.py
+++ b/src/topupopt/data/finance/invest.py
@@ -8,6 +8,8 @@ from statistics import mean
#******************************************************************************
#******************************************************************************
+# TODO: enable swapping the polarity
+
class Investment:
"""This class is meant to enable analysis of specific investments."""
@@ -174,6 +176,7 @@ class Investment:
cash_flow: float or int,
start_period: int,
longevity: int = None):
+ """Adds a sequence of cash flows to the analysis."""
if type(longevity) == type(None):
diff --git a/src/topupopt/problems/esipp/network.py b/src/topupopt/problems/esipp/network.py
index cfbb71d..e42b934 100644
--- a/src/topupopt/problems/esipp/network.py
+++ b/src/topupopt/problems/esipp/network.py
@@ -37,31 +37,20 @@ class Arcs:
validate: bool): # one number
# initialise
-
self.name = name
-
self.efficiency = efficiency
-
self.efficiency_reverse = efficiency_reverse
-
self.static_loss = static_loss
-
self.capacity = capacity
-
self.minimum_cost = minimum_cost
-
self.specific_capacity_cost = specific_capacity_cost
-
self.capacity_is_instantaneous = capacity_is_instantaneous
# results (for simulation or post-optimisation analysis)
-
self.options_selected = [False for element in self.capacity]
# validate
-
if validate:
-
Arcs.validate(self)
# *************************************************************************
@@ -148,6 +137,28 @@ class Arcs:
# *************************************************************************
# *************************************************************************
+ def has_constant_efficiency(self) -> bool:
+ """Returns True if the arc has a constant efficiency."""
+
+ if self.has_proportional_losses():
+ # proportional losses
+ if self.is_isotropic():
+ # is isotropic
+ if len(set(self.efficiency.values())) == 1:
+ # is isotropic and has constant eta = True
+ return True
+ else:
+ # is isotropic but does not have constant eta = False
+ return False
+ else: # is not isotropic
+ return False # not isotropic = efficiency can change
+ else:
+ # no proportional losses: always equal to 1 = has constant eta
+ return True
+
+ # *************************************************************************
+ # *************************************************************************
+
def has_proportional_losses(self) -> bool:
"""Returns False if the efficiency is always one, otherwise True."""
diff --git a/src/topupopt/problems/esipp/problem.py b/src/topupopt/problems/esipp/problem.py
index 9f965d1..600195c 100644
--- a/src/topupopt/problems/esipp/problem.py
+++ b/src/topupopt/problems/esipp/problem.py
@@ -29,6 +29,9 @@ from .resource import ResourcePrice
# *****************************************************************************
# *****************************************************************************
+# TODO: allow users to define how fixed components in the objective function
+# are handled (using a variable equal to one or by excluding them altogether)
+
class InfrastructurePlanningProblem(EnergySystem):
"""A class for optimisation of infrastructure planning problems."""
@@ -4225,7 +4228,8 @@ def is_peak_total_problem(problem: InfrastructurePlanningProblem) -> bool:
# conditions:
# 1) maximum congestion occurs simultaneously across the network
- # - corollary: there are no dynamic behaviours in the network
+ # - corollary: dynamic behaviours do not change the peaks
+ # - corollary: arc efficiencies are constant
# - simplifying assumption: no proportional losses in the network
# 2) the time during which maximum congestion occurs can be determined
# 3) energy prices are constant in time and volume
diff --git a/tests/examples_dhn.py b/tests/examples_dhn.py
deleted file mode 100644
index b8cd2f1..0000000
--- a/tests/examples_dhn.py
+++ /dev/null
@@ -1,1564 +0,0 @@
-# imports
-
-# standard
-
-import numbers
-
-import math
-
-import random as rand
-
-from statistics import mean
-
-# local, external
-
-import numpy as np
-
-import networkx as nx
-
-# local, internal
-
-# topupopt
-
-from src.topupopt.problems.esipp.problem import InfrastructurePlanningProblem
-
-from src.topupopt.problems.esipp.network import Network, Arcs
-
-from src.topupopt.problems.esipp.resource import ResourcePrice
-
-# import src.topupopt.problems.esipp as tuo
-
-import src.topupopt.problems.esipp.utils as utils
-
-import src.topupopt.data.dhn.utils as dhn_utils
-
-# from src.topupopt.data.dhn.network import DistrictHeatingTrenchSupplyReturnPipes
-
-from src.topupopt.data.dhn.network import PipeTrench
-
-# topupheat
-
-from topupheat.pipes.classes import StandardisedPipe, InsulatedPipe
-
-# from topupheat.pipes.trenches import SupplyReturnPipeTrenchWithIdenticalPipes
-
-from topupheat.pipes.twin import StandardisedTwinPipe, InsulatedTwinPipe
-
-# import topupheat.pipes.fic as fic
-
-import topupheat.pipes.trenches as trenches
-
-#******************************************************************************
-#******************************************************************************
-
-def examples(fluid_db,
- single_pipe_db,
- twin_pipe_db,
- solver: str,
- solver_options: dict = None,
- seed_number: int = None):
-
- #**************************************************************************
- #**************************************************************************
-
- # termination criteria
-
- solver_timelimit = 60
-
- solver_abs_mip_gap = 0.001
-
- solver_rel_mip_gal = 0.01
-
- if type(solver_options) == dict:
-
- solver_options.update({
- 'time_limit':solver_timelimit,
- 'relative_mip_gap':solver_rel_mip_gal,
- 'absolute_mip_gap':solver_abs_mip_gap
- })
-
- else:
-
- solver_options = {
- 'time_limit':solver_timelimit,
- 'relative_mip_gap':solver_rel_mip_gal,
- 'absolute_mip_gap':solver_abs_mip_gap
- }
-
- #**************************************************************************
-
- # tests
-
- validate_pipe_cost_functions()
-
- example_trench_manual(fluid_db)
-
- dht, trench_tech = example_trench_viadb(single_pipe_db, fluid_db)
-
- validate_minimum_trench_dimensions(single_pipe_db, twin_pipe_db)
-
- example_pipe_trench_objects(fluid_db, single_pipe_db, twin_pipe_db)
-
- #**************************************************************************
-
- # no sos, regular time intervals
-
- example_generic_problem_dhn(
- single_pipe_db=single_pipe_db,
- twin_pipe_db=twin_pipe_db,
- fluid_db=fluid_db,
- solver=solver,
- solver_options=solver_options,
- use_sos_arcs=False,
- sos_weight_key=InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_NONE,
- seed_number=None,
- perform_analysis=False,
- plot_results=False, # True,
- print_solver_output=False,
- irregular_time_intervals=False
- )
-
- # sos, cost as weight, regular time intervals
-
- example_generic_problem_dhn(
- single_pipe_db=single_pipe_db,
- twin_pipe_db=twin_pipe_db,
- fluid_db=fluid_db,
- solver=solver,
- solver_options=solver_options,
- use_sos_arcs=True,
- sos_weight_key=InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_COST,
- seed_number=seed_number,
- perform_analysis=False,
- plot_results=False,
- print_solver_output=False,
- irregular_time_intervals=False
- )
-
- # sos, capacity as weight, regular time intervals
-
- example_generic_problem_dhn(
- single_pipe_db=single_pipe_db,
- twin_pipe_db=twin_pipe_db,
- fluid_db=fluid_db,
- solver=solver,
- solver_options=solver_options,
- use_sos_arcs=True,
- sos_weight_key=InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_CAP,
- seed_number=seed_number,
- perform_analysis=False,
- plot_results=False,
- print_solver_output=False,
- irregular_time_intervals=False
- )
-
- # sos, specific minimum cost as weight, irregular time intervals
-
- example_generic_problem_dhn(
- single_pipe_db=single_pipe_db,
- twin_pipe_db=twin_pipe_db,
- fluid_db=fluid_db,
- solver=solver,
- solver_options=solver_options,
- use_sos_arcs=True,
- sos_weight_key=InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_SPEC_COST,
- seed_number=seed_number,
- perform_analysis=False,
- plot_results=False,
- print_solver_output=False,
- irregular_time_intervals=True
- )
-
- #**************************************************************************
- #**************************************************************************
-
-#******************************************************************************
-#******************************************************************************
-
-def example_trench_viadb(pipedb, fluiddb):
-
- # what does it do?
- # >> Creates a distric heating trench object with multiple options
-
- # seed number
-
- seed_number = 249
-
- rand.seed(seed_number)
-
- # number of intervals
-
- number_intervals = 3
-
- time_interval_duration = [rand.random() for k in range(number_intervals)]
-
- # network
-
- dhn_supply_temperature = 100+273.15 # K
-
- dhn_return_temperature = 60+273.15 # K
-
- dhn_max_specific_pressure_loss = 100 # Pa
-
- # trench
-
- trench_pipe_depth = 3
-
- trench_pipe_distance = 2 #
-
- trench_ground_thermal_conductivity = 1.8 # 0.9-2.7
-
- trench_ground_surface_temperature = [
- 7.8+273.15 for i in range(number_intervals)] # K
-
- trench_ground_air_heat_transfer_coefficient = 14.6 # W/m2
-
- # pipe details
-
- pipe_length = 1000
-
- pipe_relative_roughness = 1e-3
-
- pipe_abs_roughness = 1e-4
-
- reference_load = (50*2*pipe_length)*10
-
- #**************************************************************************
- #**************************************************************************
-
- # list_pipe_tuples = [(20,1),
- # (25,1),
- # (32,1),
- # (40,1),
- # (50,1),
- # (65,1),
- # (80,1),
- # (100,1),
- # (125,1),
- # (150,1),
- # (200,1),
- # (250,1),
- # (300,1),
- # (350,1),
- # (400,1),
- # (450,1),
- # (500,1),
- # (600,1),
- # (700,1),
- # (800,1),
- # (900,1),
- # (1000,1)]
-
- list_pipe_tuples = [pipe_tuple for pipe_tuple in pipedb.pipe_tuples]
-
- list_pipe_tuples = list_pipe_tuples[0:10]
-
- #**************************************************************************
- #**************************************************************************
-
- list_dhts = []
-
- list_pipes = []
-
- # list of pipe objects
-
- list_pipes = [
- StandardisedPipe(length=pipe_length,
- sp=0,
- pipe_tuple=pipe_tuple,
- e_rel=pipe_relative_roughness,
- # e_eff=pipe_abs_roughness,
- db=pipedb)
- for pipe_tuple in list_pipe_tuples]
-
- #**************************************************************************
- #**************************************************************************
-
- # trenches
-
- trench_tech = trenches.SupplyReturnPipeTrenchWithIdenticalPipes(
- pipes=list_pipes,
- fluid_database=fluiddb,
- ground_thermal_conductivity=trench_ground_thermal_conductivity,
- ground_air_heat_transfer_coefficient=trench_ground_air_heat_transfer_coefficient,
- pipe_center_depth=trench_pipe_depth,
- pipe_center_distance=trench_pipe_distance,
- supply_temperature=dhn_supply_temperature,
- return_temperature=dhn_return_temperature,
- max_specific_pressure_loss=dhn_max_specific_pressure_loss,
- time_interval_duration=time_interval_duration,
- surroundings_temperature=trench_ground_surface_temperature)
-
- #**************************************************************************
- #**************************************************************************
-
- # 2) create trench object
-
- dht = PipeTrench(
- name='a trench',
- trenches={0: trench_tech},
- length=pipe_length,
- use_proportional_losses=True,
- use_static_losses=False
- )
-
- #**************************************************************************
- #**************************************************************************
-
- return dht, trench_tech
-
- #**************************************************************************
- #**************************************************************************
-
-#******************************************************************************
-#******************************************************************************
-
-def example_trench_manual(fluiddb):
-
- # what does it do?
- # >> creates a district heating trench object and
-
- number_intervals = 3
-
- time_interval_duration = [rand.random() for k in range(number_intervals)]
-
- #**************************************************************************
- #**************************************************************************
-
- # network
-
- dhn_supply_temperature = 100+273.15 # K
-
- dhn_return_temperature = 60+273.15 # K
-
- dhn_max_specific_pressure_loss = 100 # Pa
-
- #**************************************************************************
- #**************************************************************************
-
- # trench
-
- trench_pipe_depth = 1.28
-
- trench_pipe_distance = 0.73
-
- trench_ground_thermal_conductivity = 0.5
-
- trench_ground_surface_temperature = [
- (7.8+273.15)+rand.random() for k in range(number_intervals)] # K
-
- trench_ground_air_heat_transfer_coefficient = 14.6 # W/m2
-
- # pipe details
-
- pipe_length = 1000
-
- pipe_k = 55
-
- pipe_e_eff = 1e-5
-
- pipe_d_int = 0.15
-
- pipe_d_ext = 0.155
-
- pipe_d_ins = 0.2
-
- pipe_d_cas = 0.205
-
- pipe_k_ins = 0.05
-
- pipe_k_cas = 0.5
-
- # pipe object
-
- pipe = StandardisedPipe(length=pipe_length,
- k=pipe_k,
- e_eff=pipe_e_eff,
- d_int=pipe_d_int,
- d_ext=pipe_d_ext,
- d_ins=pipe_d_ins,
- d_cas=pipe_d_cas,
- k_ins=pipe_k_ins,
- k_cas=pipe_k_cas,
- sp=0,
- p_rated=25)
-
- #**************************************************************************
- #**************************************************************************
-
- # trenches
-
- trench_tech = trenches.SupplyReturnPipeTrenchWithIdenticalPipes(
- pipes=[pipe],
- fluid_database=fluiddb,
- ground_thermal_conductivity=trench_ground_thermal_conductivity,
- ground_air_heat_transfer_coefficient=trench_ground_air_heat_transfer_coefficient,
- pipe_center_depth=trench_pipe_depth,
- pipe_center_distance=trench_pipe_distance,
- supply_temperature=dhn_supply_temperature,
- return_temperature=dhn_return_temperature,
- max_specific_pressure_loss=dhn_max_specific_pressure_loss,
- time_interval_duration=time_interval_duration,
- surroundings_temperature=trench_ground_surface_temperature)
-
- #**************************************************************************
- #**************************************************************************
-
- # 2) create trench object
-
- dht = PipeTrench(
- name='a trench',
- trenches={0: trench_tech},
- length=pipe_length,
- use_proportional_losses=True,
- use_static_losses=False
- )
-
- #**************************************************************************
- #**************************************************************************
-
- return dht, trench_tech
-
- #**************************************************************************
- #**************************************************************************
-
-#******************************************************************************
-#******************************************************************************
-#******************************************************************************
-#******************************************************************************
-
-def example_generic_problem_dhn(single_pipe_db,
- twin_pipe_db,
- fluid_db,
- solver: str = 'glpk',
- solver_options: dict = None,
- use_sos_arcs: bool = False,
- sos_weight_key: str = 'cost',
- seed_number: int = None,
- perform_analysis: bool = False,
- plot_results: bool = False,
- print_solver_output: bool = False,
- irregular_time_intervals: bool = False):
-
- number_periods = 2
-
- number_intraperiod_time_intervals = 5
-
- discount_rates = tuple([0.035 for p in range(number_periods)])
-
- planning_horizon = [365*24*3600 for p in range(number_periods)] # intra-period, of course
-
- if irregular_time_intervals:
-
- # TODO: adjust demand/supply levels
-
- time_step_max_relative_variation = 0.25
-
- intraperiod_time_interval_duration = [
- (planning_horizon[0]/number_intraperiod_time_intervals)*
- (1+(k/(number_intraperiod_time_intervals-1)-0.5)*
- time_step_max_relative_variation)
- for k in range(number_intraperiod_time_intervals)]
-
- else:
-
- intraperiod_time_interval_duration = [
- planning_horizon[0]/number_intraperiod_time_intervals
- for k in range(number_intraperiod_time_intervals)]
-
- # time weights
-
- # average time interval duration
-
- average_time_interval_duration = round(
- mean(
- intraperiod_time_interval_duration
- )
- )
-
- # 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
-
- # create problem object
-
- ipp = InfrastructurePlanningProblem(
- name='problem',
- discount_rates={0: discount_rates},
- reporting_periods={0: tuple(i for i in range(number_periods))},
- time_intervals={
- 0: tuple(dt for dt in intraperiod_time_interval_duration)
- },
- time_weights=time_weights
- )
-
- # add networks and systems
-
- ipp = create_generic_networks_dhn(single_pipe_db,
- twin_pipe_db,
- fluid_db,
- ipp,
- use_sos_arcs,
- sos_weight_key,
- seed_number)
-
- # set up the use of sos, if necessary
-
- 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():
-
- continue # TODO: reach this instruction
-
- ipp.use_sos1_for_arc_selection(
- network_key,
- arc_key,
- use_real_variables_if_possible=False,
- sos1_weight_method=sos_weight_key)
-
- # instantiate
-
- ipp.instantiate()
-
- # optimise
-
- ipp.optimise(solver_name=solver,
- solver_options=solver_options,
- output_options={},
- print_solver_output=print_solver_output)
-
- # run tests
-
- utils.run_mvesipp_analysis(ipp,
- ipp.instance,
- analyse_problem=perform_analysis,
- analyse_results=perform_analysis)
-
- # plot the network
-
- # summarise dhn results
-
- for key, network in ipp.networks.items():
-
- data_dict = dhn_utils.summarise_network_by_arc_technology(network)
-
- dhn_utils.print_summary_arc_technologies(data_dict)
-
- #**************************************************************************
- #**************************************************************************
-
- # print results
-
- utils.plot_networks(ipp, filepath=None, filename_radical=None)
-
- #**************************************************************************
- #**************************************************************************
-
- for key, network in ipp.networks.items():
-
- # set network crs
-
- network.graph['crs'] = "EPSG:4326"
-
- # add random network positions
-
- for node_key in network.nodes():
-
- node_dict = dict(network.nodes[node_key])
-
- node_dict['x'] = rand.random()
- node_dict['y'] = rand.random()
-
- network.modify_network_node(node_key,
- node_data=node_dict)
-
- dhn_utils.plot_network_layout(network,
- include_basemap=False)
-
- #**************************************************************************
- #**************************************************************************
-
- # return something
-
- return True
-
-#******************************************************************************
-#******************************************************************************
-#******************************************************************************
-#******************************************************************************
-
-def create_generic_networks_dhn(single_pipe_db,
- twin_pipe_db,
- fluid_db,
- ipp: InfrastructurePlanningProblem,
- use_sos_arcs: bool = False,
- sos_weight_key: str = 'cost',
- seed_number: int = None):
-
- #**************************************************************************
- #**************************************************************************
-
- # how to create a random dhn problem?
- # 1) get the raw data
- # 2) get a basic network layout with the distances
- # 3) add the extra nodes
- # 4) add extra arcs
-
- #**************************************************************************
- #**************************************************************************
-
- # set the seed number
-
- if seed_number == None:
-
- seed_number = rand.randint(1,int(1e5))
-
- print('Seed number: '+str(seed_number))
-
- # initialise random number generators
-
- rand.seed(a=seed_number)
-
- np.random.seed(seed=seed_number)
-
- #**************************************************************************
- #**************************************************************************
-
- # prices
-
- #**************************************************************************
-
- # 1 kWh = 3600*1000 J
- # X kWh = 1000*1000*1000 J
- # 1 MWh = 3600*1000*1000 J
- # X MWh = 1000*1000*1000 J
-
- kWh_DIV_MWh = 1000
-
- # GJ_DIV_kWh = 1e9/(3600*1e3)
-
- GJ_DIV_MWh = 1000/3600
-
- MWh_DIV_J = 1/(3600*1000*1000)
-
- EUR_DIV_DKK = 1/7.5
-
- relative_price_variation = 0.1
-
- #**************************************************************************
-
- # nominal gas price
-
- # 20.8287 €/GJ = 20.8287/GJ_2_kWh @ DK, 2020, household consumers
-
- # 4.6664 €/GJ = 4.6664/GJ_2_kWh @ DK, 2020, non-household consumers
-
- # Source: https://ec.europa.eu/eurostat/databrowser/view/ten00118/default/table?lang=en
-
- nominal_gas_price = 20.8287*GJ_DIV_MWh # €/MWh
-
- # nominal electricity price
-
- # 0.2833 €/kWh @ DK, 2020, household consumers
-
- # 0.0540 €/kWh @ DK, 2020, non-household consumers
-
- # Source: https://ec.europa.eu/eurostat/databrowser/view/ten00117/default/table?lang=en
-
- nominal_electricity_price = 0.2833*kWh_DIV_MWh # €/MWh
-
- # nominal heat price
-
- # 700 DKK/MWh
-
- # Source: VEKS annual report 2020: https://www.veks.dk/en/publications
-
- nominal_heat_price = 700*EUR_DIV_DKK # €/MWh
-
- # price lists
-
- # imp_gas_prices = ResourcePrice(
- # prices=[nominal_gas_price*
- # (1-relative_price_variation*rand.random())
- # for k in range(ipp.number_intraperiod_time_intervals)],
- # volumes=None
- # )
-
- # imp_elec_prices = ResourcePrice(
- # prices=[nominal_electricity_price*
- # (1-relative_price_variation*rand.random())
- # for k in range(ipp.number_intraperiod_time_intervals)],
- # volumes=None
- # )
-
- # imp_dh_prices = ResourcePrice(
- # prices=[nominal_heat_price*
- # (1-relative_price_variation*rand.random())
- # for k in range(ipp.number_intraperiod_time_intervals)],
- # volumes=None
- # )
-
- imp_gas_prices = [nominal_gas_price*
- (1-relative_price_variation*rand.random())
- for k in range(ipp.number_time_intervals[0])]
-
- imp_elec_prices = [nominal_electricity_price*
- (1-relative_price_variation*rand.random())
- for k in range(ipp.number_time_intervals[0])]
-
- imp_dh_prices = [
- ResourcePrice(
- nominal_heat_price*
- (1-relative_price_variation*rand.random())
- )
- for k in range(ipp.number_time_intervals[0])
- ]
-
- #**************************************************************************
- #**************************************************************************
-
- # heat demand
-
- dh_number_demand_nodes = 5
-
- dh_annual_demand_mwh = 11606 # MWh
-
- dh_demand_node_areas = [rand.random()
- for node in range(dh_number_demand_nodes)]
-
- dh_total_area = sum(dh_demand_node_areas)
-
- dh_demand_node_annual_demand_mwh = [
- dh_annual_demand_mwh*area/dh_total_area
- for area in dh_demand_node_areas]
-
- dh_minimum_relative_demand = [
- 0.2
- for node in range(dh_number_demand_nodes)]
-
- dh_demand_node_profiles = [
- [dh_demand_node_annual_demand_mwh[node]*
- dh_minimum_relative_demand[node]/
- ipp.number_time_intervals[0]
- +
- dh_demand_node_annual_demand_mwh[node]*
- (1-dh_minimum_relative_demand[node])*
- (1+np.sin(2*np.pi*k/ipp.number_time_intervals[0]))/
- ipp.number_time_intervals[0]
- for k in range(ipp.number_time_intervals[0])]
- for node in range(dh_number_demand_nodes)
- ]
-
- #**************************************************************************
-
- # static gas demand
-
- gas_demand = [0
- for k in range(ipp.number_time_intervals[0])]
-
- #**************************************************************************
-
- # static electrical demand
-
- elec_demand = [0
- for k in range(ipp.number_time_intervals[0])]
-
- #**************************************************************************
- #**************************************************************************
-
- # transport technologies
-
- #**************************************************************************
-
- # water pipes
-
- list_single_pipe_tuples = [pipe_tuple
- for pipe_tuple in single_pipe_db.pipe_tuples]
-
- list_twin_pipe_tuples = [pipe_tuple
- for pipe_tuple in twin_pipe_db.pipe_tuples]
-
- list_single_pipes = [
- StandardisedPipe(#length=pipe_length,
- pipe_tuple=pipe_tuple,
- #e_rel=pipe_relative_roughness,
- db=single_pipe_db)
- for i, pipe_tuple in enumerate(list_single_pipe_tuples)
- ]
-
- list_twin_pipes = [
- StandardisedTwinPipe(#length=pipe_length,
- pipe_tuple=pipe_tuple,
- #e_rel=pipe_relative_roughness,
- db=twin_pipe_db)
- for i, pipe_tuple in enumerate(list_twin_pipe_tuples)
- ]
-
- # trenches
-
- trench_ground_thermal_conductivity = 1.5 # W/mK
- trench_ground_air_heat_transfer_coefficient = 10 # W/m²K
- trench_pipe_depth = 1.5 #m ; should be adjusted for each pipe size
- trench_pipe_distance = 1.5 # m ; should be adjusted for each pipe size
- dhn_supply_temperature = 273.15+70 # Celsius
- dhn_return_temperature = 273.15+40 # Celsius
- dhn_max_specific_pressure_loss = 100 # Pa/m
- trench_ground_surface_temperature = [
- 10 for k in range(ipp.number_time_intervals[0])] # ºC
-
- trench_tech = trenches.SupplyReturnPipeTrenchWithIdenticalPipes(
- pipes=list_single_pipes,
- fluid_database=fluid_db,
- ground_thermal_conductivity=trench_ground_thermal_conductivity,
- ground_air_heat_transfer_coefficient=trench_ground_air_heat_transfer_coefficient,
- pipe_center_depth=trench_pipe_depth,
- pipe_center_distance=trench_pipe_distance,
- supply_temperature=dhn_supply_temperature,
- return_temperature=dhn_return_temperature,
- max_specific_pressure_loss=dhn_max_specific_pressure_loss,
- time_interval_duration=ipp.time_intervals[0],
- surroundings_temperature=trench_ground_surface_temperature)
-
- #**************************************************************************
-
- # gas pipes
-
- #**************************************************************************
-
- # electrical lines
-
- #**************************************************************************
- #**************************************************************************
-
- # street network
-
- number_street_network_nodes = 10
-
- G_dh = nx.binomial_tree(round(math.log2(number_street_network_nodes)),
- create_using=nx.MultiDiGraph)
-
- number_street_network_nodes = G_dh.number_of_nodes()
-
- # create copies
-
- G_elec = G_dh.copy()
-
- G_gas = G_dh.copy()
-
- # create Network objects
-
- G_dh = Network(incoming_graph_data=G_dh)
-
- G_elec = Network(incoming_graph_data=G_elec)
-
- G_gas = Network(incoming_graph_data=G_gas)
-
- #**************************************************************************
- #**************************************************************************
-
- # district heating grid
-
- #**************************************************************************
-
- # update nodes
-
- for node in G_dh.nodes():
-
- G_dh.modify_network_node(node,
- {
- G_dh.KEY_NODE_TYPE:G_dh.KEY_NODE_TYPE_WAY,
- }
- )
-
- #**************************************************************************
-
- # update arcs
-
- for arc in G_dh.edges(keys=True):
-
- trench_length = 500*rand.random() # m
-
- # new_arc_tech = DistrictHeatingTrenchSupplyReturnPipes(
- # name='fake trench',
- # efficiency=None,
- # capacity=None,
- # minimum_cost=None,
- # specific_capacity_cost=None,
- # capacity_is_instantaneous=False,
- # dht=trench_tech,
- # length=trench_length,
- # capacity_unit_conversion_factor=(
- # MWh_DIV_J*ipp.time_intervals[0][0])
- # )
-
- # new_arc_tech = DistrictHeatingTrenchSupplyReturnPipes(
- # name='fake trench',
- # trench=trench_tech,
- # length=trench_length,
- # capacity_unit_conversion_factor=(
- # MWh_DIV_J*ipp.time_intervals[0][0]
- # )
- # )
-
- new_arc_tech = PipeTrench(
- name='fake trench',
- trenches={0: trench_tech},
- length=trench_length,
- use_proportional_losses=True,
- use_static_losses=False,
- capacity_unit_conversion_factor=(
- MWh_DIV_J*ipp.time_intervals[0][0]
- )
- )
-
- G_dh.modify_network_arc(
- arc[0],
- arc[1],
- arc[2],
- data_dict={
- G_dh.KEY_ARC_TECH: new_arc_tech,
- G_dh.KEY_ARC_UND: False})
-
- #**************************************************************************
-
- # add import nodes
-
- G_dh.add_import_node(
- node_key=G_dh.number_of_nodes(),
- prices={
- (q,p,k): res_pri
- for q in range(ipp.number_assessments)
- for p in range(ipp.number_reporting_periods[q])
- for k, res_pri in enumerate(imp_dh_prices)
-
- } # key: (assessment, period, interval)
- )
-
- #**************************************************************************
-
- # add export nodes: there are no export nodes
-
- #**************************************************************************
-
- # add demand nodes
-
- for node in range(dh_number_demand_nodes):
-
- # G_dh.add_network_node(node_key=G_dh.number_of_nodes(),
- # import_prices=[],
- # export_prices=[],
- # base_flow=dh_demand_node_profiles[node]
- # )
-
- G_dh.add_source_sink_node(
- node_key=G_dh.number_of_nodes(),
- #base_flow=dh_demand_node_profiles[node]
- base_flow={
- (0, k): dh_demand_node_profiles[node][k]
- for k in range(ipp.number_time_intervals[0])
- }
- )
-
- #**************************************************************************
-
- # add arcs to ensure feasibility
-
- # for each demand node, outline a path starting from an import node
-
- demand_paths = []
-
- for node in range(dh_number_demand_nodes):
-
- # identify the initial and final nodes (may need update)
-
- initial_node = number_street_network_nodes # single import node, static
-
- final_node = number_street_network_nodes+1+node # may require updates
-
- # create a random path starting in an import node
-
- new_path = [rand.randint(0, G_dh.number_of_nodes()-1)
- for k in range(G_dh.number_of_nodes()-2)]
-
- # adjust the path to ensure it goes from A to B without loops
-
- # remove loops
-
- new_path_dummy = tuple(new_path)
-
- for node in new_path_dummy:
-
- while new_path.count(node) > 1:
-
- new_path.remove(node)
-
- # initial node
-
- try:
-
- start_pos = new_path.index(initial_node)
-
- except ValueError:
-
- # initial_node is not in the path, insert it
-
- new_path.insert(0, initial_node)
-
- start_pos = 0
-
- # final node
-
- try:
-
- end_pos = new_path.index(final_node)
-
- if start_pos > end_pos:
-
- raise ValueError
-
- except ValueError:
-
- # final_node is not in new_path, append it
-
- new_path.append(final_node)
-
- end_pos = len(new_path)-1 # zero based indexing
-
- # trim the path
-
- new_path = new_path[start_pos:end_pos+1]
-
- # add path to the list of paths
-
- demand_paths.append(new_path)
-
- #**************************************************************************
-
- # for each path, create the necessary arcs
-
- for path in demand_paths:
-
- # for each consecutive pair of nodes on the path, ensure it has arcs
-
- for path_node_pair in range(len(path)-1):
-
- node_A = path[path_node_pair]
-
- node_B = path[path_node_pair+1]
-
- if G_dh.has_edge(node_A, node_B):
-
- continue
-
- else:
-
- trench_length = 500*rand.random() # m
-
- # new_arc_tech = DistrictHeatingTrenchSupplyReturnPipes(
- # name='fake trench',
- # efficiency=None,
- # capacity=None,
- # minimum_cost=None,
- # specific_capacity_cost=None,
- # capacity_is_instantaneous=False,
- # dht=trench_tech,
- # length=trench_length,
- # capacity_unit_conversion_factor=(
- # MWh_DIV_J*
- # ipp.intraperiod_time_interval_duration[0])
- # )
-
- # new_arc_tech = DistrictHeatingTrenchSupplyReturnPipes(
- # name='fake trench',
- # trench=trench_tech,
- # length=trench_length,
- # capacity_unit_conversion_factor=(
- # MWh_DIV_J*
- # ipp.intraperiod_time_interval_duration[0])
- # )
-
- new_arc_tech = PipeTrench(
- name='fake trench',
- trenches={0: trench_tech},
- length=trench_length,
- use_proportional_losses=True,
- use_static_losses=False,
- capacity_unit_conversion_factor=(
- MWh_DIV_J*ipp.time_intervals[0][0]
- )
- )
-
- G_dh.add_directed_arc(
- node_A,
- node_B,
- arcs=new_arc_tech)
-
- #**************************************************************************
- #**************************************************************************
-
- # identify import and export nodes
-
- G_dh.identify_node_types()
-
- # add network to mves object
-
- ipp.add_network(network_key='dh',
- network=G_dh)
-
- #**************************************************************************
- #**************************************************************************
-
- return ipp
-
- #**************************************************************************
- #**************************************************************************
-
-#******************************************************************************
-#******************************************************************************
-
-def validate_pipe_cost_functions():
-
- # data from Roder et al. (2021): http://doi.org/10.5278/ijsepm.6248
-
- list_pipe_dn = [25,32,40,50,63,75,90,110,125,160]
-
- list_true_specific_cost = [
- 91.046,
- 107.079,
- 127.179,
- 152.063,
- 187.712,
- 222.316,
- 268.43,
- 332.77,
- 384.748,
- 511.279]
-
- rel_tol = 0.01
-
- for index, pipe_dn in enumerate(list_pipe_dn):
-
- np.testing.assert_allclose(
- dhn_utils.specific_pipe_cost_function_thermos(pipe_dn/1000),
- list_true_specific_cost[index],
- rtol=rel_tol)
-
-#******************************************************************************
-#******************************************************************************
-
-# trench dimensions
-
-def validate_minimum_trench_dimensions(single_pipe_db,
- twin_pipe_db):
-
- # water pipes
-
- list_single_pipe_tuples = [pipe_tuple
- for pipe_tuple in single_pipe_db.pipe_tuples]
-
- list_twin_pipe_tuples = [pipe_tuple
- for pipe_tuple in twin_pipe_db.pipe_tuples]
-
- list_single_pipes = [
- StandardisedPipe(#length=pipe_length,
- pipe_tuple=pipe_tuple,
- #e_rel=pipe_relative_roughness,
- db=single_pipe_db)
- for i, pipe_tuple in enumerate(list_single_pipe_tuples)
- ]
-
- list_twin_pipes = [
- StandardisedTwinPipe(#length=pipe_length,
- pipe_tuple=pipe_tuple,
- #e_rel=pipe_relative_roughness,
- db=twin_pipe_db)
- for i, pipe_tuple in enumerate(list_twin_pipe_tuples)
- ]
-
- # add small pipes to evaluate all possible cases
-
- small_pipe = InsulatedPipe(length=6,
- k=100,
- k_ins=0.1,
- k_cas=5,
- e_eff=1e-4,
- d_int=5e-3,
- d_ext=6e-3,
- d_ins=7e-3,
- d_cas=8e-3)
-
- list_single_pipes.append(small_pipe)
-
- medium_pipe = InsulatedPipe(length=6,
- k=100,
- k_ins=0.1,
- k_cas=5,
- e_eff=1e-3,
- d_int=45e-3,
- d_ext=50e-3,
- d_ins=80e-3,
- d_cas=88e-3)
-
- list_single_pipes.append(medium_pipe)
-
- small_twin_pipe = InsulatedTwinPipe(length=6,
- k=100,
- k_ins=0.1,
- k_cas=5,
- e_eff=1e-4,
- d_int=2e-3,
- d_ext=2.5e-3,
- d_ins=8.0e-3,
- d_cas=8.5e-3,
- pipe_center_distance=3.5e-3)
-
- list_twin_pipes.append(small_twin_pipe)
-
- medium_twin_pipe = InsulatedTwinPipe(length=6,
- k=100,
- k_ins=0.1,
- k_cas=5,
- e_eff=1e-3,
- d_int=2e-3,
- d_ext=2.5e-3,
- d_ins=75e-3,
- d_cas=80e-3,
- pipe_center_distance=3.5e-3)
-
- list_twin_pipes.append(medium_twin_pipe)
-
- #**************************************************************************
-
- list_pipes = []
- list_pipes.extend(list_single_pipes)
- list_pipes.extend(list_twin_pipes)
-
- # given list_pipes
-
- for pipe in list_pipes:
-
- if isinstance(pipe, InsulatedTwinPipe):
-
- pipe.ValidateInsulatedTwinPipe()
-
- pipe_depth = dhn_utils.recommended_minimum_pipe_center_depth(pipe)
-
- assert pipe_depth > 0
-
- assert pipe_depth > pipe.d_cas/2
-
- elif isinstance(pipe, InsulatedPipe):
-
- pipe.ValidateInsulatedPipe()
-
- pipe_depth = dhn_utils.recommended_minimum_pipe_center_depth(pipe)
-
- assert pipe_depth > 0
-
- assert pipe_depth > pipe.d_cas/2
-
- pipe_distance = dhn_utils.recommended_minimum_interpipe_distance(
- pipe)
-
- assert pipe_distance > 0
-
- assert pipe_distance > pipe.d_cas
-
- # else:
-
- # raise NotImplementedError
-
-#******************************************************************************
-#******************************************************************************
-
-# test pipe trench objects
-
-def example_pipe_trench_objects(fluid_db,
- single_pipe_db,
- twin_pipe_db):
-
- #**************************************************************************
- #**************************************************************************
-
- # water pipes
-
- list_single_pipe_tuples = [pipe_tuple
- for pipe_tuple in single_pipe_db.pipe_tuples]
-
- list_twin_pipe_tuples = [pipe_tuple
- for pipe_tuple in twin_pipe_db.pipe_tuples]
-
- list_single_pipes = [
- StandardisedPipe(pipe_tuple=pipe_tuple,
- db=single_pipe_db)
- for i, pipe_tuple in enumerate(list_single_pipe_tuples)
- ]
-
- list_twin_pipes = [
- StandardisedTwinPipe(pipe_tuple=pipe_tuple,
- db=twin_pipe_db)
- for i, pipe_tuple in enumerate(list_twin_pipe_tuples)
- ]
-
- #**************************************************************************
-
- # what does it do?
- # >> Creates a distric heating trench object with multiple options
-
- # seed number
-
- seed_number = 249
-
- rand.seed(seed_number)
-
- # number of intervals
-
- number_intervals = 3
-
- time_interval_duration = [rand.random() for k in range(number_intervals)]
-
- # network
-
- dhn_supply_temperature = 100+273.15 # K
-
- dhn_return_temperature = 60+273.15 # K
-
- dhn_max_specific_pressure_loss = 100 # Pa
-
- # trench
-
- trench_pipe_depth = 3
-
- trench_pipe_distance = 2 #
-
- trench_ground_thermal_conductivity = 1.8 # 0.9-2.7
-
- trench_ground_surface_temperature = [
- 7.8+273.15 for i in range(number_intervals)] # K
-
- trench_ground_air_heat_transfer_coefficient = 14.6 # W/m2
-
- # pipe details
-
- pipe_length = 1000
-
- pipe_relative_roughness = 1e-3
-
- #**************************************************************************
- #**************************************************************************
-
- # single pipe trenches
-
- trench_tech = trenches.SupplyReturnPipeTrenchWithIdenticalPipes(
- pipes=list_single_pipes,
- fluid_database=fluid_db,
- ground_thermal_conductivity=trench_ground_thermal_conductivity,
- ground_air_heat_transfer_coefficient=trench_ground_air_heat_transfer_coefficient,
- pipe_center_depth=trench_pipe_depth,
- pipe_center_distance=trench_pipe_distance,
- supply_temperature=dhn_supply_temperature,
- return_temperature=dhn_return_temperature,
- max_specific_pressure_loss=dhn_max_specific_pressure_loss,
- time_interval_duration=time_interval_duration,
- surroundings_temperature=trench_ground_surface_temperature)
-
- # single pipe, no external cost, no offset
-
- pipe_trench_obj = PipeTrench(name='hello',
- trenches={0: trench_tech},
- length=pipe_length,
- use_proportional_losses=True,
- use_static_losses=True,
- minimum_cost=None,
- minimum_cost_offset=None,
- validate=True)
-
- original_min_cost = tuple(pipe_trench_obj.minimum_cost)
-
- # single pipe, no external cost, offset
-
- pipe_trench_obj = PipeTrench(name='hello',
- trenches={0: trench_tech},
- length=pipe_length,
- use_proportional_losses=True,
- use_static_losses=True,
- minimum_cost=None,
- minimum_cost_offset=tuple(
- rand.random()
- for pipe in list_single_pipes
- ),
- validate=True)
-
- for orig_cost, new_cost in zip(original_min_cost,
- pipe_trench_obj.minimum_cost):
-
- assert orig_cost <= new_cost
-
- # single pipe, external cost, no offset
-
- external_cost = tuple(0.2+min_cost for min_cost in original_min_cost)
-
- pipe_trench_obj = PipeTrench(name='hello',
- trenches={0: trench_tech},
- length=pipe_length,
- use_proportional_losses=True,
- use_static_losses=True,
- minimum_cost=external_cost,
- minimum_cost_offset=None,
- validate=True)
-
- assert external_cost == pipe_trench_obj.minimum_cost
-
- # single pipe, external cost, offset
-
- error_triggered = False
- try:
- pipe_trench_obj = PipeTrench(name='hello',
- trenches={0: trench_tech},
- length=pipe_length,
- use_proportional_losses=True,
- use_static_losses=True,
- minimum_cost=original_min_cost,
- minimum_cost_offset=external_cost,
- validate=True)
- except TypeError:
- error_triggered = True
- assert error_triggered
-
- # use list as minimum cost offset
-
- error_triggered = False
- try:
- pipe_trench_obj = PipeTrench(name='hello',
- trenches={0: trench_tech},
- length=pipe_length,
- use_proportional_losses=True,
- use_static_losses=True,
- minimum_cost=None,
- minimum_cost_offset=list(
- rand.random()
- for pipe in list_single_pipes
- ),
- validate=True)
- except TypeError:
- error_triggered = True
- assert error_triggered
-
- #**************************************************************************
- #**************************************************************************
-
- # twin pipe trenches
-
- trench_tech = trenches.SupplyReturnPipeTrenchWithIdenticalPipes(
- pipes=list_twin_pipes,
- fluid_database=fluid_db,
- ground_thermal_conductivity=trench_ground_thermal_conductivity,
- ground_air_heat_transfer_coefficient=trench_ground_air_heat_transfer_coefficient,
- pipe_center_depth=trench_pipe_depth,
- pipe_center_distance=trench_pipe_distance,
- supply_temperature=dhn_supply_temperature,
- return_temperature=dhn_return_temperature,
- max_specific_pressure_loss=dhn_max_specific_pressure_loss,
- time_interval_duration=time_interval_duration,
- surroundings_temperature=trench_ground_surface_temperature)
-
- # single pipe, no external cost, no offset
-
- pipe_trench_obj = PipeTrench(name='hello',
- trenches={0: trench_tech},
- length=pipe_length,
- use_proportional_losses=True,
- use_static_losses=True,
- minimum_cost=None,
- minimum_cost_offset=None,
- validate=True)
-
- original_min_cost = tuple(pipe_trench_obj.minimum_cost)
-
- # single pipe, no external cost, offset
-
- pipe_trench_obj = PipeTrench(name='hello',
- trenches={0: trench_tech},
- length=pipe_length,
- use_proportional_losses=True,
- use_static_losses=True,
- minimum_cost=None,
- minimum_cost_offset=tuple(
- rand.random()
- for pipe in list_twin_pipes
- ),
- validate=True)
-
- for orig_cost, new_cost in zip(original_min_cost,
- pipe_trench_obj.minimum_cost):
-
- assert orig_cost <= new_cost
-
- # single pipe, external cost, no offset
-
- external_cost = tuple(0.2+min_cost for min_cost in original_min_cost)
-
- pipe_trench_obj = PipeTrench(name='hello',
- trenches={0: trench_tech},
- length=pipe_length,
- use_proportional_losses=True,
- use_static_losses=True,
- minimum_cost=external_cost,
- minimum_cost_offset=None,
- validate=True)
-
- assert external_cost == pipe_trench_obj.minimum_cost
-
- # single pipe, external cost, offset
-
- error_triggered = False
- try:
- pipe_trench_obj = PipeTrench(name='hello',
- trenches={0: trench_tech},
- length=pipe_length,
- use_proportional_losses=True,
- use_static_losses=True,
- minimum_cost=original_min_cost,
- minimum_cost_offset=external_cost,
- validate=True)
- except TypeError:
- error_triggered = True
- assert error_triggered
-
- # use list as minimum cost offset
-
- error_triggered = False
- try:
- pipe_trench_obj = PipeTrench(name='hello',
- trenches={0: trench_tech},
- length=pipe_length,
- use_proportional_losses=True,
- use_static_losses=True,
- minimum_cost=None,
- minimum_cost_offset=list(
- rand.random()
- for pipe in list_twin_pipes
- ),
- validate=True)
- except TypeError:
- error_triggered = True
- assert error_triggered
-
- #**************************************************************************
- #**************************************************************************
-
-#******************************************************************************
-#******************************************************************************
\ No newline at end of file
diff --git a/tests/examples_dynsys.py b/tests/examples_dynsys.py
deleted file mode 100644
index dd9c404..0000000
--- a/tests/examples_dynsys.py
+++ /dev/null
@@ -1,2366 +0,0 @@
-# imports
-
-# standard
-
-import random
-
-# local, external
-
-import numpy as np
-
-from scipy.signal import StateSpace, lsim, dlsim
-
-# local, internal
-
-import src.topupopt.problems.esipp.dynsys as dynsys
-
-#******************************************************************************
-#******************************************************************************
-
-def examples():
-
- #**************************************************************************
- #**************************************************************************
-
- seed_number = random.randint(1,int(1e5))
-
- print_outputs = True
-
- # with states and outputs
-
- # test multi-ODE, multi-output dynamic systems while integrating outputs
-
- examples_dynsys_multiode_multiout(True, print_outputs, seed_number)
-
- # test multi-ODE, multi-output dynamic systems without integrating outputs
-
- examples_dynsys_multiode_multiout(False, print_outputs, seed_number)
-
- # test single ODE, multi-output dynamic systems while integrating outputs
-
- examples_dynsys_singleode_multiout(True, print_outputs, seed_number)
-
- # test single ODE, multi-output dynamic systems without integrating outputs
-
- examples_dynsys_singleode_multiout(False, print_outputs, seed_number)
-
- # test multi-ODE, single-output dynamic systems while integrating outputs
-
- examples_dynsys_multiode_multiout(True, print_outputs, seed_number, 1)
-
- # test multi-ODE, single-output dynamic systems without integrating outputs
-
- examples_dynsys_multiode_multiout(False, print_outputs, seed_number, 1)
-
- # test single-ODE, single-output dynamic systems while integrating outputs
-
- examples_dynsys_singleode_multiout(True, print_outputs, seed_number, 1)
-
- # test single-ODE, single-output dynamic systems without integrating outputs
-
- examples_dynsys_singleode_multiout(False, print_outputs, seed_number, 1)
-
- #**************************************************************************
-
- # outputless
-
- # test single-ODE, outputless dynamic systems while integrating outputs
-
- examples_dynsys_singleode_multiout(True, print_outputs, seed_number, 0)
-
- # test multi-ODE, outputless dynamic systems while integrating outputs
-
- examples_dynsys_multiode_multiout(True, print_outputs, seed_number, 0)
-
- # test single-ODE, outputless dynamic systems without integrating outputs
-
- examples_dynsys_singleode_multiout(False, print_outputs, seed_number, 0)
-
- # test multi-ODE, outputless dynamic systems without integrating outputs
-
- examples_dynsys_multiode_multiout(False, print_outputs, seed_number, 0)
-
- # outputless system via dynsys subclass
-
- example_outputless_system_object()
-
- #**************************************************************************
-
- # stateless
-
- # test stateless, single-output dynamic systems while integrating outputs
-
- examples_dynsys_stateless_multiout(True, print_outputs, seed_number, 1)
-
- # test stateless, multi-output dynamic systems without integrating outputs
-
- examples_dynsys_stateless_multiout(False, print_outputs, seed_number, 2)
-
- # stateless system via dynsys subclass
-
- example_stateless_system_object(True)
- example_stateless_system_object(False)
-
- #**************************************************************************
- #**************************************************************************
-
- # trigger errors
-
- # test stateless, outputless dynamic systems while integrating outputs
-
- number_errors = 0
-
- try:
- examples_dynsys_stateless_multiout(True, False, seed_number, 0)
- except Exception:
- number_errors += 1
-
- assert number_errors == 1
-
- # test stateless, outputless dynamic systems without integrating outputs
-
- number_errors = 0
-
- try:
- examples_dynsys_stateless_multiout(False, False, seed_number, 0)
- except Exception:
- number_errors += 1
-
- assert number_errors == 1
-
- # test negative time duration
-
- example_incorrect_time_step_durations()
-
- # test unrecognised matrix formats
-
- example_unrecognised_matrix_formats()
-
- # different matrix sizes for the same problem, all other things being equal
-
- example_varying_matrix_sizes(True)
- example_varying_matrix_sizes(False)
-
- # test multiple A matrices and multiple non-matching time intervals
-
- example_nonmatching_time_steps_and_matrices()
-
- # test non-square A matrices
-
- example_nonsquare_A_matrices()
-
- # test incompatible A and B matrices (different number of rows)
-
- example_incompatible_AB_matrices()
-
- # test incompatible C and D matrices (different number of rows)
-
- example_incompatible_CD_matrices()
-
- # test incompatible A and C matrices (different number of columns)
-
- example_incompatible_AC_matrices()
-
- # test incompatible B and D matrices (different number of columns)
-
- example_incompatible_BD_matrices()
-
- # trigger incorrect input signal format error when simulating
-
- example_single_time_step_model_incorrect_inputs()
-
- # TODO: test only some matrices as being time invariant
-
- #**************************************************************************
- #**************************************************************************
-
-#******************************************************************************
-#******************************************************************************
-
-def example_incorrect_time_step_durations(integrate_outputs: bool = True):
-
- #**************************************************************************
-
- # negative time step duration
-
- # regular time steps (as an int)
-
- time_step_durations = [-1]
-
- # get the model
-
- Aw, min_rel_heat = get_stateless_model_data()
-
- _, _, _, D = stateless_model(Aw, min_rel_heat)
-
- number_errors = 0
- try:
- _ = dynsys.StatelessSystem(
- time_interval_durations=time_step_durations[0],
- D=D,
- integrate_outputs=integrate_outputs)
- except ValueError:
- number_errors += 1
-
- assert number_errors == 1
-
- #**************************************************************************
-
- # no time step duration (empty list)
-
- time_step_durations = []
-
- # get the model
-
- Aw, min_rel_heat = get_stateless_model_data()
-
- _, _, _, D = stateless_model(Aw, min_rel_heat)
-
- number_errors = 0
- try:
- _ = dynsys.StatelessSystem(
- time_interval_durations=time_step_durations,
- D=D,
- integrate_outputs=integrate_outputs)
- except ValueError:
- number_errors += 1
-
- assert number_errors == 1
-
- #**************************************************************************
-
- # time step duration list with negative or zero time step durations
-
- time_step_durations = [-1, 3, 0]
-
- # get the model
-
- Aw, min_rel_heat = get_stateless_model_data()
-
- _, _, _, D = stateless_model(Aw, min_rel_heat)
-
- number_errors = 0
- try:
- _ = dynsys.StatelessSystem(
- time_interval_durations=time_step_durations,
- D=D,
- integrate_outputs=integrate_outputs)
- except ValueError:
- number_errors += 1
-
- assert number_errors == 1
-
- #**************************************************************************
-
- # time step duration list with non-numeric types
-
- time_step_durations = [None, 3, 3]
-
- # get the model
-
- Aw, min_rel_heat = get_stateless_model_data()
-
- _, _, _, D = stateless_model(Aw, min_rel_heat)
-
- number_errors = 0
- try:
- _ = dynsys.StatelessSystem(
- time_interval_durations=time_step_durations,
- D=D,
- integrate_outputs=integrate_outputs)
- except TypeError:
- number_errors += 1
-
- assert number_errors == 1
-
-#******************************************************************************
-#******************************************************************************
-
-def example_single_time_step_model_incorrect_inputs(integrate_y: bool = True):
-
- # regular time steps (as an int)
-
- time_step_durations = [1]
-
- # define a sequence of time steps
-
- t = np.array([sum(time_step_durations[0:i])
- for i in range(len(time_step_durations)+1)])
-
- # define the inputs
-
- list_inputs = [-5, 50, 0.1, 1, 5] # extra input signal to force error
-
- U = np.array([[input_i for dt in t] # extra time step to force special case
- for input_i in list_inputs])
-
- # U = np.array([[input_i for _ in range(len(time_step_durations))]
- # for input_i in list_inputs])
-
- # get the model
-
- Ci, Ch, Ria, Rih, Aw, min_rel_heat, x0 = get_multi_ode_model_data()
-
- A, B, C, D = two_node_model(Ci, Ch, Ria, Rih, Aw, min_rel_heat)
-
- ds = dynsys.DynamicSystem(
- time_interval_durations=time_step_durations,
- A=A,
- B=B,
- C=C,
- D=D,
- integrate_outputs=integrate_y)
-
- # define the initial conditions
-
- number_errors = 0
-
- try:
- X, Y = ds.simulate(U, x0)
- except ValueError:
- number_errors += 1
-
- assert number_errors == 1
-
-#******************************************************************************
-#******************************************************************************
-
-def example_outputless_system_object():
-
- # regular time steps
-
- time_step_durations = [1, 1, 1, 1]
-
- # define a sequence of time steps
-
- t = np.array([sum(time_step_durations[0:i])
- for i in range(len(time_step_durations)+1)])
-
- # define the inputs
-
- list_inputs = [-5, 50, 0.1, 1]
-
- U = np.array([[input_i for _ in range(len(time_step_durations))]
- for input_i in list_inputs])
-
- # get the model
-
- Ci, Ch, Ria, Rih, Aw, min_rel_heat, x0 = get_multi_ode_model_data()
-
- A, B, _, _ = two_node_model(Ci, Ch, Ria, Rih, Aw, min_rel_heat)
-
- ds = dynsys.OutputlessSystem(
- time_interval_durations=time_step_durations,
- A=A,
- B=B)
-
- # define the initial conditions
-
- X, Y = ds.simulate(U, x0)
-
- assert Y == None
-
- assert isinstance(X, np.ndarray)
-
- #**************************************************************************
-
- # regular time steps (as an int)
-
- time_step_durations = [1]
-
- # define a sequence of time steps
-
- t = np.array([sum(time_step_durations[0:i])
- for i in range(len(time_step_durations)+1)])
-
- # define the inputs
-
- list_inputs = [-5, 50, 0.1, 1]
-
- U = np.array([[input_i for _ in range(len(time_step_durations))]
- for input_i in list_inputs])
-
- # get the model
-
- Ci, Ch, Ria, Rih, Aw, min_rel_heat, x0 = get_multi_ode_model_data()
-
- A, B, _, _ = two_node_model(Ci, Ch, Ria, Rih, Aw, min_rel_heat)
-
- ds = dynsys.OutputlessSystem(
- time_interval_durations=time_step_durations,
- A=A,
- B=B)
-
- # define the initial conditions
-
- X, Y = ds.simulate(U, x0)
-
- assert Y == None
-
- assert isinstance(X, np.ndarray)
-
- #**************************************************************************
-
- # irregular time steps
-
- time_step_durations = [1, 1.5, 0.5, 1]
-
- # define a sequence of time steps
-
- t = np.array([sum(time_step_durations[0:i])
- for i in range(len(time_step_durations)+1)])
-
- # define the inputs
-
- list_inputs = [-5, 50, 0.1, 1]
-
- U = np.array([[input_i for _ in range(len(time_step_durations))]
- for input_i in list_inputs])
-
- # get the model
-
- Ci, Ch, Ria, Rih, Aw, min_rel_heat, x0 = get_multi_ode_model_data()
-
- A, B, _, _ = two_node_model(Ci, Ch, Ria, Rih, Aw, min_rel_heat)
-
- ds = dynsys.OutputlessSystem(
- time_interval_durations=time_step_durations,
- A=A,
- B=B)
-
- # define the initial conditions
-
- X, Y = ds.simulate(U, x0)
-
- assert Y == None
-
- assert isinstance(X, np.ndarray)
-
-#******************************************************************************
-#******************************************************************************
-
-def example_stateless_system_object(integrate_outputs: bool = True):
-
- # regular time steps
-
- time_step_durations = [1, 1, 1, 1]
-
- # define a sequence of time steps
-
- t = np.array([sum(time_step_durations[0:i])
- for i in range(len(time_step_durations)+1)])
-
- # define the inputs
-
- list_inputs = [-5, 50, 0.1, 1]
-
- U = np.array([[input_i for _ in range(len(time_step_durations))]
- for input_i in list_inputs])
-
- # get the model
-
- Aw, min_rel_heat = get_stateless_model_data()
-
- _, _, _, D = stateless_model(Aw, min_rel_heat)
-
- ds = dynsys.StatelessSystem(
- time_interval_durations=time_step_durations,
- D=D,
- integrate_outputs=integrate_outputs)
-
- # define the initial conditions
-
- X, Y = ds.simulate(U)
-
- assert X == None
-
- assert isinstance(Y, np.ndarray)
-
- if integrate_outputs:
-
- assert Y.shape == (2,len(t)-1) # two outputs
-
- else:
-
- assert Y.shape == (2,len(t)) # two outputs
-
- #**************************************************************************
-
- # regular time steps (as an int)
-
- time_step_durations = [1]
-
- # define a sequence of time steps
-
- t = np.array([sum(time_step_durations[0:i])
- for i in range(len(time_step_durations)+1)])
-
- # define the inputs
-
- list_inputs = [-5, 50, 0.1, 1]
-
- U = np.array([[input_i for _ in range(len(time_step_durations))]
- for input_i in list_inputs])
-
- # get the model
-
- Aw, min_rel_heat = get_stateless_model_data()
-
- _, _, _, D = stateless_model(Aw, min_rel_heat)
-
- ds = dynsys.StatelessSystem(
- time_interval_durations=time_step_durations[0],
- D=D,
- integrate_outputs=integrate_outputs)
-
- # define the initial conditions
-
- X, Y = ds.simulate(U)
-
- # X, Y = ds.simulate(U[:,:-1]) # to avoid having the second
-
- assert X == None
-
- assert isinstance(Y, np.ndarray)
-
- if integrate_outputs:
-
- assert Y.shape == (2,len(t)-1) # two outputs
-
- else:
-
- assert Y.shape == (2,len(t)) # two outputs
-
- #**************************************************************************
-
- # irregular time steps
-
- time_step_durations = [1, 1.5, 0.5, 1]
-
- # define a sequence of time steps
-
- t = np.array([sum(time_step_durations[0:i])
- for i in range(len(time_step_durations)+1)])
-
- # define the inputs
-
- list_inputs = [-5, 50, 0.1, 1]
-
- U = np.array([[input_i for _ in range(len(time_step_durations))]
- for input_i in list_inputs])
-
- # get the model
-
- Aw, min_rel_heat = get_stateless_model_data()
-
- _, _, _, D = stateless_model(Aw, min_rel_heat)
-
- ds = dynsys.StatelessSystem(
- time_interval_durations=time_step_durations,
- D=D,
- integrate_outputs=integrate_outputs)
-
- # define the initial conditions
-
- X, Y = ds.simulate(U)
-
- assert X == None
-
- assert isinstance(Y, np.ndarray)
-
- if integrate_outputs:
-
- assert Y.shape == (2,len(t)-1) # two outputs
-
- else:
-
- assert Y.shape == (2,len(t)) # two outputs
-
-#******************************************************************************
-#******************************************************************************
-
-def example_incompatible_BD_matrices(integrate_y: bool = True):
-
- # B and D must have the same number of columns
-
- time_step_durations = [1, 1, 1]
-
- A = np.array([[4, 9], [1, 3]])
-
- B = np.array([[5, 6, 1], [7, 2, 1]])
-
- C = np.array([[4, 6]])
-
- D = np.array([[3, 6]])
-
- # test
-
- number_errors = 0
-
- try:
- _ = dynsys.DynamicSystem(
- time_interval_durations=time_step_durations,
- A=A,
- B=B,
- C=C,
- D=D,
- integrate_outputs=integrate_y)
- except ValueError:
- number_errors += 1
-
- assert number_errors == 1
-
-#******************************************************************************
-#******************************************************************************
-
-def example_incompatible_AC_matrices(integrate_y: bool = True):
-
- # A and C must have the same number of columns
-
- time_step_durations = [1, 1, 1]
-
- A = np.array([[4, 9], [1, 3]])
-
- B = np.array([[5, 6, 1], [7, 2, 1]])
-
- C = np.array([[4]])
-
- D = np.array([[3, 6, 3]])
-
- # test
-
- number_errors = 0
-
- try:
- _ = dynsys.DynamicSystem(
- time_interval_durations=time_step_durations,
- A=A,
- B=B,
- C=C,
- D=D,
- integrate_outputs=integrate_y)
- except ValueError:
- number_errors += 1
-
- assert number_errors == 1
-
-#******************************************************************************
-#******************************************************************************
-
-def example_incompatible_CD_matrices(integrate_y: bool = True):
-
- # C and D must have the same number of rows
-
- time_step_durations = [1, 1, 8]
-
- A = np.array([[4, 9], [1, 3]])
-
- B = np.array([[5, 6, 1], [7, 2, 1]])
-
- C = np.array([[4, 9], [1, 3]])
-
- D = np.array([[3, 6, 3]])
-
- # test
-
- number_errors = 0
-
- try:
- _ = dynsys.DynamicSystem(
- time_interval_durations=time_step_durations,
- A=A,
- B=B,
- C=C,
- D=D,
- integrate_outputs=integrate_y)
- except ValueError:
- number_errors += 1
-
- assert number_errors == 1
-
-#******************************************************************************
-#******************************************************************************
-
-def example_incompatible_AB_matrices(integrate_y: bool = True):
-
- # A and B must have the same number of rows
-
- time_step_durations = [1, 1, 8]
-
- A = np.array([[4, 9],[1, 3]])
-
- B = np.array([[3, 6, 3]])
-
- C = np.array([[2, 4]])
-
- D = np.array([[5, 6, 1]])
-
- # test A
-
- number_errors = 0
-
- try:
- _ = dynsys.DynamicSystem(
- time_interval_durations=time_step_durations,
- A=A,
- B=B,
- C=C,
- D=D,
- integrate_outputs=integrate_y)
- except ValueError:
- number_errors += 1
-
- assert number_errors == 1
-
-#******************************************************************************
-#******************************************************************************
-
-def example_nonsquare_A_matrices(integrate_y: bool = True):
-
- # Single non-square A matrix
-
- time_step_durations = [1, 1, 8]
-
- A = np.array([[4, 9]])
-
- B = np.array([[3]])
-
- C = np.array([[2]])
-
- D = np.array([[5]])
-
- # test A
-
- number_errors = 0
-
- try:
- _ = dynsys.DynamicSystem(
- time_interval_durations=time_step_durations,
- A=A,
- B=B,
- C=C,
- D=D,
- integrate_outputs=integrate_y)
- except ValueError:
- number_errors += 1
-
- assert number_errors == 1
-
- # # multiple matrices, at least one of which is non-square
-
- # A = [np.array([[4]]), np.array([[1, 9]])]
-
- # B = [np.array([[3]]), np.array([[4]])]
-
- # C = [np.array([[2]]), np.array([[2]])]
-
- # D = [np.array([[8]]), np.array([[7]])]
-
- # # test A
-
- # number_errors = 0
-
- # try:
- # _ = dynsys.DynamicSystem(
- # time_interval_durations=time_step_durations,
- # A=A,
- # B=B,
- # C=C,
- # D=D,
- # integrate_outputs=integrate_y)
- # except ValueError:
- # number_errors += 1
-
- # assert number_errors == 1
-
-#******************************************************************************
-#******************************************************************************
-
-def example_nonmatching_time_steps_and_matrices(integrate_y: bool = True):
-
- # multiple matrices and time intervals, but more time steps than matrices
-
- time_step_durations = [1, 1, 8]
-
- A = [np.array([[4]]),np.array([[1]])]
-
- B = [np.array([[3]]),np.array([[4]]),np.array([[4]])]
-
- C = [np.array([[2]]),np.array([[2]]),np.array([[2]])]
-
- D = [np.array([[8]]),np.array([[7]]),np.array([[7]])]
-
- # test A
-
- number_errors = 0
-
- try:
- _ = dynsys.DynamicSystem(
- time_interval_durations=time_step_durations,
- A=A,
- B=B,
- C=C,
- D=D,
- integrate_outputs=integrate_y)
- except ValueError:
- number_errors += 1
-
- assert number_errors == 1
-
- # test B
-
- number_errors = 0
-
- try:
- _ = dynsys.DynamicSystem(
- time_interval_durations=time_step_durations,
- A=B,
- B=A,
- C=C,
- D=D,
- integrate_outputs=integrate_y)
- except ValueError:
- number_errors += 1
-
- assert number_errors == 1
-
- # test C
-
- number_errors = 0
-
- try:
- _ = dynsys.DynamicSystem(
- time_interval_durations=time_step_durations,
- A=C,
- B=B,
- C=A,
- D=D,
- integrate_outputs=integrate_y)
- except ValueError:
- number_errors += 1
-
- assert number_errors == 1
-
- # test D
-
- number_errors = 0
-
- try:
- _ = dynsys.DynamicSystem(
- time_interval_durations=time_step_durations,
- A=D,
- B=B,
- C=C,
- D=A,
- integrate_outputs=integrate_y)
- except ValueError:
- number_errors += 1
-
- assert number_errors == 1
-
- # multiple matrices and time intervals, but more matrices than time steps
-
-#******************************************************************************
-#******************************************************************************
-
-def example_unrecognised_matrix_formats(integrate_outputs: bool = True):
-
- # single matrix: matrix as lists of lists
-
- time_step_durations = 8
-
- A = 3
-
- B = np.array([[3]])
-
- C = np.array([[2]])
-
- D = np.array([[8]])
-
- # test A
-
- number_errors = 0
-
- try:
- _ = dynsys.DynamicSystem(
- time_interval_durations=time_step_durations,
- A=A,
- B=B,
- C=C,
- D=D,
- integrate_outputs=integrate_outputs)
- except TypeError:
- number_errors += 1
-
- assert number_errors == 1
-
- # test B
-
- number_errors = 0
-
- try:
- _ = dynsys.DynamicSystem(
- time_interval_durations=time_step_durations,
- A=B,
- B=A,
- C=C,
- D=D,
- integrate_outputs=integrate_outputs)
- except TypeError:
- number_errors += 1
-
- assert number_errors == 1
-
- # test C
-
- number_errors = 0
-
- try:
- _ = dynsys.DynamicSystem(
- time_interval_durations=time_step_durations,
- A=C,
- B=B,
- C=A,
- D=D,
- integrate_outputs=integrate_outputs)
- except TypeError:
- number_errors += 1
-
- assert number_errors == 1
-
- # test D
-
- number_errors = 0
-
- try:
- _ = dynsys.DynamicSystem(
- time_interval_durations=time_step_durations,
- A=D,
- B=B,
- C=C,
- D=A,
- integrate_outputs=integrate_outputs)
- except TypeError:
- number_errors += 1
-
- assert number_errors == 1
-
- #**************************************************************************
-
- # multiple matrices: matrix as lists of lists
-
- time_step_durations = [1, 1]
-
- A = [[[3]], [[3]]]
-
- B = [np.array([[3]]),np.array([[4]])]
-
- C = [np.array([[2]]),np.array([[2]])]
-
- D = [np.array([[8]]),np.array([[7]])]
-
- # test A
-
- number_errors = 0
-
- try:
- _ = dynsys.DynamicSystem(
- time_interval_durations=time_step_durations,
- A=A,
- B=B,
- C=C,
- D=D,
- integrate_outputs=integrate_outputs)
- except TypeError:
- number_errors += 1
-
- assert number_errors == 1
-
- # test B
-
- number_errors = 0
-
- try:
- _ = dynsys.DynamicSystem(
- time_interval_durations=time_step_durations,
- A=B,
- B=A,
- C=C,
- D=D,
- integrate_outputs=integrate_outputs)
- except TypeError:
- number_errors += 1
-
- assert number_errors == 1
-
- # test C
-
- number_errors = 0
-
- try:
- _ = dynsys.DynamicSystem(
- time_interval_durations=time_step_durations,
- A=C,
- B=B,
- C=A,
- D=D,
- integrate_outputs=integrate_outputs)
- except TypeError:
- number_errors += 1
-
- assert number_errors == 1
-
- # test D
-
- number_errors = 0
-
- try:
- _ = dynsys.DynamicSystem(
- time_interval_durations=time_step_durations,
- A=D,
- B=B,
- C=C,
- D=A,
- integrate_outputs=integrate_outputs)
- except TypeError:
- number_errors += 1
-
- assert number_errors == 1
-
- #**************************************************************************
-
-#******************************************************************************
-#******************************************************************************
-
-def example_varying_matrix_sizes(integrate_outputs: bool):
-
- time_step_durations = [1, 1]
-
- A1 = np.array([[1]])
-
- A2 = np.array([[1,3],[4,7]])
-
- A = [A1, A2]
-
- B = [np.array([[3]]),np.array([[4]])]
-
- C = [np.array([[2]]),np.array([[2]])]
-
- D = [np.array([[8]]),np.array([[7]])]
-
- # test A
-
- number_errors = 0
-
- try:
- _ = dynsys.DynamicSystem(
- time_interval_durations=time_step_durations,
- A=A,
- B=B,
- C=C,
- D=D,
- integrate_outputs=integrate_outputs)
- except ValueError:
- number_errors += 1
-
- assert number_errors == 1
-
- # test B
-
- number_errors = 0
-
- try:
- _ = dynsys.DynamicSystem(
- time_interval_durations=time_step_durations,
- A=B,
- B=A,
- C=C,
- D=D,
- integrate_outputs=integrate_outputs)
- except ValueError:
- number_errors += 1
-
- assert number_errors == 1
-
- # test C
-
- number_errors = 0
-
- try:
- _ = dynsys.DynamicSystem(
- time_interval_durations=time_step_durations,
- A=C,
- B=B,
- C=A,
- D=D,
- integrate_outputs=integrate_outputs)
- except ValueError:
- number_errors += 1
-
- assert number_errors == 1
-
- # test D
-
- number_errors = 0
-
- try:
- _ = dynsys.DynamicSystem(
- time_interval_durations=time_step_durations,
- A=D,
- B=B,
- C=C,
- D=A,
- integrate_outputs=integrate_outputs)
- except ValueError:
- number_errors += 1
-
- assert number_errors == 1
-
-#******************************************************************************
-#******************************************************************************
-
-def examples_dynsys_stateless_multiout(integrate_outputs: bool,
- print_output: bool = True,
- seed_number: int = None,
- number_outputs: int = 2):
-
- if number_outputs == 1:
-
- print('******************************************************************')
- print('******************************************************************')
- print('stateless, single output')
-
- else:
-
- print('******************************************************************')
- print('******************************************************************')
- print('stateless, multi output')
-
- #**************************************************************************
-
- # time steps
-
- list_regular_time_steps = [1, 1, 1, 1]
-
- list_irregular_time_steps = [1, 1.5, 0.5, 1]
-
- assert len(list_regular_time_steps) == len(list_irregular_time_steps)
-
- # define the inputs
-
- list_inputs = [-5, 50, 0.1, 1]
-
- #**************************************************************************
-
- # generate time varying problem
-
- list_A_matrices = []
-
- list_B_matrices = []
-
- list_C_matrices = []
-
- list_D_matrices = []
-
- for dt in list_irregular_time_steps:
-
- # data
-
- Aw, min_rel_heat = get_stateless_model_data(
- relative_amplitude_variation=0.1)
-
- # matrices
-
- (A_matrix,
- B_matrix,
- C_matrix,
- D_matrix) = stateless_model(Aw, min_rel_heat)
-
- if number_outputs == 0:
-
- D_matrix = None
-
- elif number_outputs != None:
-
- D_matrix = D_matrix[0:number_outputs,:]
-
- list_A_matrices.append(A_matrix)
- list_B_matrices.append(B_matrix)
- list_C_matrices.append(C_matrix)
- list_D_matrices.append(D_matrix)
-
- if number_outputs == 0:
-
- list_D_matrices = None
-
- #**************************************************************************
-
- # generate time invariant problem
-
- # data
-
- Aw, min_rel_heat = get_stateless_model_data()
-
- # matrices
-
- (A_matrix,
- B_matrix,
- C_matrix,
- D_matrix) = stateless_model(Aw, min_rel_heat)
-
- if number_outputs == 0:
-
- D_matrix = None
-
- elif number_outputs != None:
-
- D_matrix = D_matrix[0:number_outputs,:]
-
- #**************************************************************************
-
- # time invariant model, regular time steps
-
- # x_scipy, y_scipy = example_scipy_time_invariant_regular_steps(
- # list_regular_time_steps,
- # list_inputs,
- # integrate_outputs,
- # A_matrix,
- # B_matrix,
- # C_matrix,
- # D_matrix,
- # x0=None)
-
- x, y = example_dynsys_time_invariant(
- list_regular_time_steps,
- list_inputs,
- integrate_outputs,
- A_matrix,
- B_matrix,
- C_matrix,
- D_matrix,
- x0=None)
-
- if print_output:
-
- print('time invariant model, regular time steps')
-
- # print(x_scipy)
- print(x)
- # print(y_scipy)
- print(y)
-
- # for (x_scipy_i, x_i) in zip(x_scipy, x):
-
- # np.testing.assert_allclose(x_i, x_scipy_i)
-
- # if number_outputs != 0:
-
- # for (y_scipy_i, y_i) in zip(y_scipy, y):
-
- # np.testing.assert_allclose(y_i, y_scipy_i)
-
- #**************************************************************************
-
- # time invariant model, irregular time steps
-
- # x_scipy, y_scipy = example_scipy_time_invariant_irregular_steps(
- # list_irregular_time_steps,
- # list_inputs,
- # integrate_outputs,
- # A_matrix,
- # B_matrix,
- # C_matrix,
- # D_matrix,
- # x0=None)
-
- x, y = example_dynsys_time_invariant(
- list_irregular_time_steps,
- list_inputs,
- integrate_outputs,
- A_matrix,
- B_matrix,
- C_matrix,
- D_matrix,
- x0=None)
-
- if print_output:
-
- print('time invariant model, irregular time steps')
-
- # print(x_scipy)
- print(x)
- # print(y_scipy)
- print(y)
-
- # for (x_scipy_i, x_i) in zip(x_scipy, x):
-
- # np.testing.assert_allclose(x_i, x_scipy_i)
-
- # if number_outputs != 0:
-
- # for (y_scipy_i, y_i) in zip(y_scipy, y):
-
- # np.testing.assert_allclose(y_i, y_scipy_i)
-
- #**************************************************************************
-
-#******************************************************************************
-#******************************************************************************
-
-def examples_dynsys_multiode_multiout(integrate_outputs: bool,
- print_output: bool = True,
- seed_number: int = None,
- number_outputs: int = 2):
-
- if number_outputs == 0:
-
- print('******************************************************************')
- print('******************************************************************')
- print('multi ODE, outputless')
-
- elif number_outputs == 1:
-
- print('******************************************************************')
- print('******************************************************************')
- print('multi ODE, single output')
-
- else:
-
- print('******************************************************************')
- print('******************************************************************')
- print('multi ODE, multi output')
-
- #**************************************************************************
-
- # time steps
-
- list_regular_time_steps = [1, 1, 1, 1]
-
- list_irregular_time_steps = [1, 1.5, 0.5, 1]
-
- assert len(list_regular_time_steps) == len(list_irregular_time_steps)
-
- # define the inputs
-
- list_inputs = [-5, 50, 0.1, 1]
-
- #**************************************************************************
-
- # generate time varying problem
-
- list_A_matrices = []
-
- list_B_matrices = []
-
- list_C_matrices = []
-
- list_D_matrices = []
-
- for dt in list_irregular_time_steps:
-
- # data
-
- Ci, Ch, Ria, Rih, Aw, min_rel_heat, x0 = get_multi_ode_model_data(
- relative_amplitude_variation=0.1)
-
- # matrices
-
- (A_matrix,
- B_matrix,
- C_matrix,
- D_matrix) = two_node_model(Ci, Ch, Ria, Rih, Aw, min_rel_heat)
-
- if number_outputs == 0:
-
- C_matrix = None
-
- D_matrix = None
-
- elif number_outputs != None:
-
- C_matrix = C_matrix[0:number_outputs,:]
-
- D_matrix = D_matrix[0:number_outputs,:]
-
- list_A_matrices.append(A_matrix)
- list_B_matrices.append(B_matrix)
- list_C_matrices.append(C_matrix)
- list_D_matrices.append(D_matrix)
-
- if number_outputs == 0:
-
- list_C_matrices = None
-
- list_D_matrices = None
-
- #**************************************************************************
-
- # generate time invariant problem
-
- # data
-
- Ci, Ch, Ria, Rih, Aw, min_rel_heat, x0 = get_multi_ode_model_data()
-
- # matrices
-
- (A_matrix,
- B_matrix,
- C_matrix,
- D_matrix) = two_node_model(Ci, Ch, Ria, Rih, Aw, min_rel_heat)
-
- if number_outputs == 0:
-
- C_matrix = None
-
- D_matrix = None
-
- elif number_outputs != None:
-
- C_matrix = C_matrix[0:number_outputs,:]
-
- D_matrix = D_matrix[0:number_outputs,:]
-
- #**************************************************************************
-
- # time invariant model, regular time steps
-
- x_scipy, y_scipy = example_scipy_time_invariant_regular_steps(
- list_regular_time_steps,
- list_inputs,
- integrate_outputs,
- A_matrix,
- B_matrix,
- C_matrix,
- D_matrix,
- x0)
-
- x, y = example_dynsys_time_invariant(
- list_regular_time_steps,
- list_inputs,
- integrate_outputs,
- A_matrix,
- B_matrix,
- C_matrix,
- D_matrix,
- x0)
-
- if print_output:
-
- print('time invariant model, regular time steps')
-
- print(x_scipy)
- print(x)
- print(y_scipy)
- print(y)
-
- for (x_scipy_i, x_i) in zip(x_scipy, x):
-
- np.testing.assert_allclose(x_i, x_scipy_i)
-
- if number_outputs != 0:
-
- for (y_scipy_i, y_i) in zip(y_scipy, y):
-
- np.testing.assert_allclose(y_i, y_scipy_i)
-
- #**************************************************************************
-
- # time invariant model, irregular time steps
-
- x_scipy, y_scipy = example_scipy_time_invariant_irregular_steps(
- list_irregular_time_steps,
- list_inputs,
- integrate_outputs,
- A_matrix,
- B_matrix,
- C_matrix,
- D_matrix,
- x0)
-
- x, y = example_dynsys_time_invariant(
- list_irregular_time_steps,
- list_inputs,
- integrate_outputs,
- A_matrix,
- B_matrix,
- C_matrix,
- D_matrix,
- x0)
-
- if print_output:
-
- print('time invariant model, irregular time steps')
-
- print(x_scipy)
- print(x)
- print(y_scipy)
- print(y)
-
- for (x_scipy_i, x_i) in zip(x_scipy, x):
-
- np.testing.assert_allclose(x_i, x_scipy_i)
-
- if number_outputs != 0:
-
- for (y_scipy_i, y_i) in zip(y_scipy, y):
-
- np.testing.assert_allclose(y_i, y_scipy_i)
-
- #**************************************************************************
-
- # time-varying model, regular time steps
-
- x_scipy, y_scipy = example_scipy_time_varying(
- list_regular_time_steps,
- list_inputs,
- integrate_outputs,
- list_A_matrices,
- list_B_matrices,
- list_C_matrices,
- list_D_matrices,
- x0)
-
- x, y = example_dynsys_time_varying(
- list_regular_time_steps,
- list_inputs,
- integrate_outputs,
- list_A_matrices,
- list_B_matrices,
- list_C_matrices,
- list_D_matrices,
- x0)
-
- if print_output:
-
- print('time-varying model, regular time steps')
-
- print(x_scipy)
- print(x)
- print(y_scipy)
- print(y)
-
- for (x_scipy_i, x_i) in zip(x_scipy, x):
-
- np.testing.assert_allclose(x_i, x_scipy_i)
-
- if number_outputs != 0:
-
- for (y_scipy_i, y_i) in zip(y_scipy, y):
-
- np.testing.assert_allclose(y_i, y_scipy_i)
-
- #**************************************************************************
-
- # time-varying model, irregular time steps
-
- x_scipy, y_scipy = example_scipy_time_varying(
- list_irregular_time_steps,
- list_inputs,
- integrate_outputs,
- list_A_matrices,
- list_B_matrices,
- list_C_matrices,
- list_D_matrices,
- x0)
-
- x, y = example_dynsys_time_varying(
- list_irregular_time_steps,
- list_inputs,
- integrate_outputs,
- list_A_matrices,
- list_B_matrices,
- list_C_matrices,
- list_D_matrices,
- x0)
-
- if print_output:
-
- print('time-varying model, irregular time steps')
-
- print(x_scipy)
- print(x)
- print(y_scipy)
- print(y)
-
- for (x_scipy_i, x_i) in zip(x_scipy, x):
-
- np.testing.assert_allclose(x_i, x_scipy_i)
-
- if number_outputs != 0:
-
- for (y_scipy_i, y_i) in zip(y_scipy, y):
-
- np.testing.assert_allclose(y_i, y_scipy_i)
-
-#******************************************************************************
-#******************************************************************************
-
-def examples_dynsys_singleode_multiout(integrate_outputs: bool,
- print_output: bool = True,
- seed_number: int = None,
- number_outputs: int = 2):
-
- if number_outputs == 0:
-
- print('******************************************************************')
- print('******************************************************************')
- print('single ODE, outputless')
-
- elif number_outputs == 1:
-
- print('******************************************************************')
- print('******************************************************************')
- print('single ODE, single output')
-
- else:
-
- print('******************************************************************')
- print('******************************************************************')
- print('single ODE, multi output')
-
- #**************************************************************************
-
- # time steps
-
- list_regular_time_steps = [1, 1, 1, 1]
-
- list_irregular_time_steps = [1, 1.5, 0.5, 1]
-
- assert len(list_regular_time_steps) == len(list_irregular_time_steps)
-
- # define the inputs
-
- list_inputs = [-5, 50, 0.1, 1]
-
- #**************************************************************************
-
- # generate time varying problem
-
- list_A_matrices = []
-
- list_B_matrices = []
-
- list_C_matrices = []
-
- list_D_matrices = []
-
- for dt in list_irregular_time_steps:
-
- # data
-
- Ci, Ria, Aw, min_rel_heat, x0 = get_single_ode_model_data(
- relative_amplitude_variation=0.1)
-
- # matrices
-
- (A_matrix,
- B_matrix,
- C_matrix,
- D_matrix) = single_node_model(Ci, Ria, Aw, min_rel_heat)
-
- if number_outputs == 0:
-
- C_matrix = None
-
- D_matrix = None
-
- elif number_outputs != None:
-
- C_matrix = C_matrix[0:number_outputs,:]
-
- D_matrix = D_matrix[0:number_outputs,:]
-
- list_A_matrices.append(A_matrix)
- list_B_matrices.append(B_matrix)
- list_C_matrices.append(C_matrix)
- list_D_matrices.append(D_matrix)
-
- if number_outputs == 0:
-
- list_C_matrices = None
-
- list_D_matrices = None
-
- #**************************************************************************
-
- # generate time invariant problem
-
- # data
-
- Ci, Ria, Aw, min_rel_heat, x0 = get_single_ode_model_data()
-
- # matrices
-
- (A_matrix,
- B_matrix,
- C_matrix,
- D_matrix) = single_node_model(Ci, Ria, Aw, min_rel_heat)
-
- if number_outputs == 0:
-
- C_matrix = None
-
- D_matrix = None
-
- elif number_outputs != None:
-
- C_matrix = C_matrix[0:number_outputs,:]
-
- D_matrix = D_matrix[0:number_outputs,:]
-
- #**************************************************************************
-
- # time invariant model, regular time steps
-
- x_scipy, y_scipy = example_scipy_time_invariant_regular_steps(
- list_regular_time_steps,
- list_inputs,
- integrate_outputs,
- A_matrix,
- B_matrix,
- C_matrix,
- D_matrix,
- x0)
-
- x, y = example_dynsys_time_invariant(
- list_regular_time_steps,
- list_inputs,
- integrate_outputs,
- A_matrix,
- B_matrix,
- C_matrix,
- D_matrix,
- x0)
-
- if print_output:
-
- print('time invariant model, regular time steps')
-
- print(x_scipy)
- print(x)
- print(y_scipy)
- print(y)
-
- for (x_scipy_i, x_i) in zip(x_scipy, x):
-
- np.testing.assert_allclose(x_i, x_scipy_i)
-
- if number_outputs != 0:
-
- for (y_scipy_i, y_i) in zip(y_scipy, y):
-
- np.testing.assert_allclose(y_i, y_scipy_i)
-
- #**************************************************************************
-
- # time invariant model, irregular time steps
-
- x_scipy, y_scipy = example_scipy_time_invariant_irregular_steps(
- list_irregular_time_steps,
- list_inputs,
- integrate_outputs,
- A_matrix,
- B_matrix,
- C_matrix,
- D_matrix,
- x0)
-
- x, y = example_dynsys_time_invariant(
- list_irregular_time_steps,
- list_inputs,
- integrate_outputs,
- A_matrix,
- B_matrix,
- C_matrix,
- D_matrix,
- x0)
-
- if print_output:
-
- print('time invariant model, irregular time steps')
-
- print(x_scipy)
- print(x)
- print(y_scipy)
- print(y)
-
- for (x_scipy_i, x_i) in zip(x_scipy, x):
-
- np.testing.assert_allclose(x_i, x_scipy_i)
-
- if number_outputs != 0:
-
- for (y_scipy_i, y_i) in zip(y_scipy, y):
-
- np.testing.assert_allclose(y_i, y_scipy_i)
-
- #**************************************************************************
-
- # time-varying model, regular time steps
-
- x_scipy, y_scipy = example_scipy_time_varying(
- list_regular_time_steps,
- list_inputs,
- integrate_outputs,
- list_A_matrices,
- list_B_matrices,
- list_C_matrices,
- list_D_matrices,
- x0)
-
- x, y = example_dynsys_time_varying(
- list_regular_time_steps,
- list_inputs,
- integrate_outputs,
- list_A_matrices,
- list_B_matrices,
- list_C_matrices,
- list_D_matrices,
- x0)
-
- if print_output:
-
- print('time-varying model, regular time steps')
-
- print(x_scipy)
- print(x)
- print(y_scipy)
- print(y)
-
- for (x_scipy_i, x_i) in zip(x_scipy, x):
-
- np.testing.assert_allclose(x_i, x_scipy_i)
-
- if number_outputs != 0:
-
- for (y_scipy_i, y_i) in zip(y_scipy, y):
-
- np.testing.assert_allclose(y_i, y_scipy_i)
-
- #**************************************************************************
-
- # time-varying model, irregular time steps
-
- x_scipy, y_scipy = example_scipy_time_varying(
- list_irregular_time_steps,
- list_inputs,
- integrate_outputs,
- list_A_matrices,
- list_B_matrices,
- list_C_matrices,
- list_D_matrices,
- x0)
-
- x, y = example_dynsys_time_varying(
- list_irregular_time_steps,
- list_inputs,
- integrate_outputs,
- list_A_matrices,
- list_B_matrices,
- list_C_matrices,
- list_D_matrices,
- x0)
-
- if print_output:
-
- print('time-varying model, irregular time steps')
-
- print(x_scipy)
- print(x)
- print(y_scipy)
- print(y)
-
- for (x_scipy_i, x_i) in zip(x_scipy, x):
-
- np.testing.assert_allclose(x_i, x_scipy_i)
-
- if number_outputs != 0:
-
- for (y_scipy_i, y_i) in zip(y_scipy, y):
-
- np.testing.assert_allclose(y_i, y_scipy_i)
-
-#******************************************************************************
-#******************************************************************************
-
-def example_scipy_time_invariant_regular_steps(
- time_step_durations: list,
- inputs_list: list,
- integrate_outputs: bool,
- A_matrix,
- B_matrix,
- C_matrix,
- D_matrix,
- x0):
-
- # define a sequence of time steps
-
- t = np.array([sum(time_step_durations[0:i])
- for i in range(len(time_step_durations)+1)])
-
- # define the inputs
-
- u = np.array([inputs_list for dt in t])
-
- if type(C_matrix) == type(None) and type(D_matrix) == type(None):
-
- ss = StateSpace(A_matrix,
- B_matrix,
- np.zeros((1,A_matrix.shape[1])),
- np.zeros((1,B_matrix.shape[1])))
-
- else:
-
- ss = StateSpace(A_matrix, B_matrix, C_matrix, D_matrix)
-
- # simulate the system response
-
- tout, yout, xout = lsim(ss, U=u, T=t, X0=x0)
-
- # convert to matrix format if there is only one dimension
-
- if len(xout.shape) == 1:
-
- xout = np.array([xout]).T
-
- if len(yout.shape) == 1:
-
- yout = np.array([yout]).T
-
- if integrate_outputs:
-
- # yout's shape should be: (number of time steps, 2)
-
- # ignore the first instant
-
- yout = yout[1:,:]
-
- yout_m, yout_n = yout.shape
-
- # multiply each output by the respective time step
-
- yout = np.array(
- [[yout[m,n]*time_step_durations[n] for n in range(yout_n)]
- for m in range(yout_m)]
- )
-
- # note: the code above should not work when C != 0
-
- return xout.T, yout.T
-
-#******************************************************************************
-#******************************************************************************
-
-def example_scipy_time_invariant_irregular_steps(
- time_step_durations: list,
- inputs_list: list,
- integrate_outputs: bool,
- A_matrix,
- B_matrix,
- C_matrix,
- D_matrix,
- x0):
-
- # define a sequence of time steps
-
- t = np.array([sum(time_step_durations[0:i])
- for i in range(len(time_step_durations)+1)])
-
- # define the inputs
-
- u = np.array([inputs_list for _ in t])
-
- # get the model
-
- if type(C_matrix) == type(None) and type(D_matrix) == type(None):
-
- ss = StateSpace(A_matrix,
- B_matrix,
- np.zeros((1,A_matrix.shape[1])),
- np.zeros((1,B_matrix.shape[1])))
-
- yout = np.zeros((1, len(time_step_durations)+1))
-
- else:
-
- ss = StateSpace(A_matrix, B_matrix, C_matrix, D_matrix)
-
- yout = np.zeros((C_matrix.shape[0], len(time_step_durations)+1))
-
- # declare output variables
-
- xout = np.zeros((A_matrix.shape[0], len(time_step_durations)+1))
-
- #yout = np.zeros((C_matrix.shape[0], len(time_step_durations)+1))
-
- xout[:,0] = x0
-
- # initial output
-
- if not integrate_outputs:
-
- # do not ignore the first instant
-
- yout[:,0] = np.dot(ss.D, u[0,:])
-
- # else:
-
- # yout[:,0] = np.array([0, 0])
-
- # for each time step
-
- for i in range(len(time_step_durations)):
-
- # time step duration
-
- dt = time_step_durations[i]
-
- ssd = ss.to_discrete(dt=dt)
-
- # simulate the system response
-
- _, ytemp, xtemp = dlsim(ssd,
- u=np.array([u[i,:],u[i+1,:]]),
- t=np.array([0,dt]),
- x0=xout[:,i])
-
- xout[:,i+1] = xtemp[1,:] #np.array([0,0])
-
- yout[:,i+1] = ytemp[1,:] #np.array([0,0])
-
- if integrate_outputs:
-
- # yout's shape should be: (2,number of time steps)
-
- # ignore the first instant
-
- yout = yout[:,1:]
-
- yout_m, yout_n = yout.shape
-
- # multiply each output by the respective time step
-
- for m in range(yout_m):
-
- for n in range(yout_n):
-
- if n == 0:
-
- continue # ignore the first time interval
-
- yout[m,n] = yout[m,n]*time_step_durations[n]
-
- return xout, yout
-
-#******************************************************************************
-#******************************************************************************
-
-def example_scipy_time_varying(time_step_durations: list,
- inputs_list: list,
- integrate_outputs: bool,
- A_matrix,
- B_matrix,
- C_matrix,
- D_matrix,
- x0):
-
- # define a sequence of time steps
-
- t = np.array([sum(time_step_durations[0:i])
- for i in range(len(time_step_durations)+1)])
-
- # define the inputs
-
- u = np.array([inputs_list for dt in t])
-
- # define the initial conditions
-
- xout = np.zeros((A_matrix[0].shape[0], len(time_step_durations)+1))
-
- xout[:,0] = x0
-
- if type(C_matrix) == type(None) and type(D_matrix) == type(None):
-
- yout = np.zeros((1, len(time_step_durations)+1))
-
- else:
-
- yout = np.zeros((C_matrix[0].shape[0], len(time_step_durations)+1))
-
-
- for i, dt in enumerate(time_step_durations):
-
- # get the model
-
- if type(C_matrix) == type(None) and type(D_matrix) == type(None):
-
- ss = StateSpace(A_matrix[i],
- B_matrix[i],
- np.zeros((1,A_matrix[i].shape[1])),
- np.zeros((1,B_matrix[i].shape[1])))
-
- else:
-
- ss = StateSpace(A_matrix[i], B_matrix[i], C_matrix[i], D_matrix[i])
-
- if i == 0 and not integrate_outputs:
-
- # compute the initial output
-
- yout[:,0] = np.dot(ss.D, u[0,:])
-
- ssd = ss.to_discrete(dt=dt)
-
- # simulate the system response
-
- if i == 0:
-
- _, ytemp, xtemp = dlsim(ssd,
- u=np.array([u[i],u[i]]),
- t=np.array([0,dt]),
- x0=x0)
-
- else:
-
- _, ytemp, xtemp = dlsim(ssd,
- u=np.array([u[i],u[i]]),
- t=np.array([0,dt]),
- x0=xtemp[1,:])
-
- # assert that the sizes match
-
- xout[:,i+1] = xtemp[1,:]
-
- yout[:,i+1] = ytemp[1,:]
-
- if integrate_outputs:
-
- # yout's shape should be: (2,number of time steps)
-
- # ignore the first instant
-
- yout = yout[:,1:]
-
- yout_m, yout_n = yout.shape
-
- # multiply each output by the respective time step
-
- for m in range(yout_m):
-
- for n in range(yout_n):
-
- if n == 0:
-
- continue # ignore the first time interval
-
- yout[m,n] = yout[m,n]*time_step_durations[n]
-
- return xout, yout
-
-#******************************************************************************
-#******************************************************************************
-
-def example_dynsys_time_invariant(time_step_durations: list,
- inputs_list: list,
- integrate_outputs: bool,
- A_matrix,
- B_matrix,
- C_matrix,
- D_matrix,
- x0):
-
- # define a sequence of time steps
-
- t = np.array([sum(time_step_durations[0:i])
- for i in range(len(time_step_durations)+1)])
-
- # define the inputs
-
- # U = np.array([[input_i for dt in t]
- # for input_i in inputs_list])
-
- U = np.array([[input_i for _ in range(len(time_step_durations))]
- for input_i in inputs_list])
-
- # get the model
-
- ds = dynsys.DynamicSystem(
- time_interval_durations=time_step_durations,
- A=A_matrix,
- B=B_matrix,
- C=C_matrix,
- D=D_matrix,
- integrate_outputs=integrate_outputs)
-
- # define the initial conditions
-
- X, Y = ds.simulate(U, x0)
-
- return X, Y
-
-#******************************************************************************
-#******************************************************************************
-
-def example_dynsys_time_varying(time_step_durations: list,
- inputs_list: list,
- integrate_outputs: bool,
- A_matrix,
- B_matrix,
- C_matrix,
- D_matrix,
- x0):
-
- # define a sequence of time steps
-
- t = np.array([sum(time_step_durations[0:i])
- for i in range(len(time_step_durations)+1)])
-
- # define the inputs
-
- # U = np.array([[input_i for dt in t]
- # for input_i in inputs_list])
-
- U = np.array([[input_i for _ in range(len(time_step_durations))]
- for input_i in inputs_list])
-
- ds = dynsys.DynamicSystem(
- time_interval_durations=time_step_durations,
- A=A_matrix,
- B=B_matrix,
- C=C_matrix,
- D=D_matrix,
- integrate_outputs=integrate_outputs)
-
- # define the initial conditions
-
- X, Y = ds.simulate(U, x0)
-
- return X, Y
-
-#******************************************************************************
-#******************************************************************************
-
-def get_stateless_model_data(relative_amplitude_variation: float = 0.0):
-
- mrh_deviation = random.random() - 0.5
-
- Aw = 6.22 # original: 6.22 m2
-
- min_rel_heat = 0.2*(1+relative_amplitude_variation*mrh_deviation)
-
- return Aw, min_rel_heat
-
-#******************************************************************************
-#******************************************************************************
-
-def get_single_ode_model_data(relative_amplitude_variation: float = 0.0):
-
- # define how the coefficients change
-
- Ria_deviation = random.random() - 0.5
-
- # define the (A, B, C and D) matrices
- # A: n*n
- # B: n*m
- # C: r*n
- # D: r*m
-
- Ci = 1.360*3600000
- Ria = (
- (1+relative_amplitude_variation*Ria_deviation)*5.31/3600000
- )
- Aw = 6.22
-
- min_rel_heat = 0.2
-
- x0 = np.array([20])
-
- return Ci, Ria, Aw, min_rel_heat, x0
-
-#******************************************************************************
-#******************************************************************************
-
-def get_multi_ode_model_data(relative_amplitude_variation: float = 0.0):
-
- # define how the coefficients change
-
- Rih_deviation = random.random() - 0.5
-
- Ria_deviation = random.random() - 0.5
-
- # define the (A, B, C and D) matrices
- # A: n*n
- # B: n*m
- # C: r*n
- # D: r*m
-
- # from Bacher and Madsen (2011): model TiTh
-
- Ci = 1.360*3600000 # original: 1.36 kWh/ºC
- Ch = 0.309*3600000 # original: 0.309 kWh/ºC
- Ria = (
- (1+relative_amplitude_variation*Ria_deviation)*5.31/3600000
- ) # original: 5.31 ºC/kWh
- Rih = (
- (1+relative_amplitude_variation*Rih_deviation)*0.639/3600000
- ) # original: 0.639 ºC/kWh
- Aw = 6.22 # original: 6.22 m2
-
- min_rel_heat = 0.2
-
- x0 = np.array([20, 20])
-
- return Ci, Ch, Ria, Rih, Aw, min_rel_heat, x0
-
-#******************************************************************************
-#******************************************************************************
-
-def stateless_model(Aw, min_rel_heat):
-
- # inputs: Ta, phi_s, phi_h above the minimum, phi_h status
- # outputs: solar irradiance, heat
-
- d = np.array([[0, Aw, 0, 0],
- [0, 0, (1-min_rel_heat), min_rel_heat]])
-
- return None, None, None, d
-
-#******************************************************************************
-#******************************************************************************
-
-def single_node_model(Ci, Ria, Aw, min_rel_heat):
-
- # states: Ti and Th
- # inputs: Ta, phi_s, phi_h above the minimum, phi_h status
- # outputs: solar irradiance, heat
-
- a = np.array([[-1/(Ria*Ci)]])
- b = np.array([[1/(Ci*Ria), Aw/Ci, (1-min_rel_heat)/Ci, min_rel_heat/Ci]])
- c = np.array([[0],[0]])
- d = np.array([[0, Aw, 0, 0],
- [0, 0, (1-min_rel_heat), min_rel_heat]])
-
- return a, b, c, d
-
-#******************************************************************************
-#******************************************************************************
-
-def two_node_model(Ci, Ch, Ria, Rih, Aw, min_rel_heat):
-
- # states: Ti and Th
- # inputs: Ta, phi_s, phi_h above the minimum, phi_h status
- # outputs: solar irradiance, heat
-
- a = np.array([[-(1/Rih+1/Ria)/Ci, 1/(Ci*Rih)],
- [1/(Ch*Rih), -1/(Ch*Rih)]])
- b = np.array([[1/(Ci*Ria), Aw/Ci, 0, 0],
- [0, 0, (1-min_rel_heat)/Ch, min_rel_heat/Ch]])
- c = np.array([[0, 0],
- [0, 0]])
- d = np.array([[0, Aw, 0, 0],
- [0, 0, (1-min_rel_heat), min_rel_heat]])
-
- return a, b, c, d
-
-#******************************************************************************
-#******************************************************************************
\ No newline at end of file
diff --git a/tests/test_all.py b/tests/test_all.py
index ca908e5..4cae3ee 100644
--- a/tests/test_all.py
+++ b/tests/test_all.py
@@ -5,9 +5,6 @@ import random
from topupheat.pipes.single import StandardisedPipeDatabase
from topupheat.common.fluids import FluidDatabase#, Fluid
-from examples_converter import examples as examples_converter
-# from examples_dhn import examples as examples_dhn
-from examples_dynsys import examples as examples_dynsys
from examples_esipp_network import examples as examples_esipp_network
from examples_esipp_problem import examples as examples_esipp_problem
from examples_esipp import examples as examples_esipp
@@ -20,9 +17,6 @@ from examples_signal import examples as examples_signal
def test_suite():
- test_examples_converter = True
- # test_examples_converter = False
-
test_examples_dynsys = True
# test_examples_dynsys = False
@@ -135,30 +129,6 @@ def test_suite():
#**************************************************************************
#**************************************************************************
- # converter
-
- if test_examples_converter:
-
- print('\'converter\': testing about to start...')
-
- examples_converter()
-
- print('\'converter\': testing complete.')
-
- #**************************************************************************
-
- # dynsys
-
- if test_examples_dynsys:
-
- print('\'dynsys\': testing about to start...')
-
- examples_dynsys()
-
- print('\'dynsys\': testing complete.')
-
- #**************************************************************************
-
# esipp-network
if test_examples_esipp_network:
diff --git a/tests/test_converter.py b/tests/test_converter.py
deleted file mode 100644
index 5850885..0000000
--- a/tests/test_converter.py
+++ /dev/null
@@ -1,292 +0,0 @@
-# imports
-
-# standard
-
-import random
-
-# local, external
-
-import numpy as np
-
-# local, internal
-
-import src.topupopt.problems.esipp.dynsys as dynsys
-
-import src.topupopt.problems.esipp.converter as cvn
-
-import src.topupopt.problems.esipp.signal as sgn
-
-# *****************************************************************************
-# *****************************************************************************
-
-class TestConverter:
-
- # test converters
- # 1) regular and irregular time steps
- # 2) time invariant and time varying models
- # 3) integrate and do not integrate outputs
- # 4) generate coefficients
-
- # test creating a stateless converter without outputs
-
- # test creating a stateless converter with 1 output
-
- # test creating a stateless converter with 2 output
-
- # test creating a converter based on a single ODE system without outputs
-
- # test creating a converter based on a single ODE system with 1 output
-
- # test creating a converter based on a single ODE system with 2 output
-
- # test creating a converter based on a multiple ODE system without outputs
-
- # test creating a converter based on a multiple ODE system with 1 output
-
- # test creating a converter based on a multiple ODE multi-output system
-
- def test_converter_regular_timesteps(self):
- time_step_durations = [1,1,1,1]
- self.example_full_converter(time_step_durations)
-
- def test_converter_irregular_timesteps(self):
- time_step_durations = [1,1.5,0.5,1]
- self.example_full_converter(time_step_durations)
-
- # *************************************************************************
- # *************************************************************************
-
- def get_stateless_model_data(relative_amplitude_variation: float = 0.0):
-
- mrh_deviation = random.random() - 0.5
-
- Aw = 6.22 # original: 6.22 m2
-
- min_rel_heat = 0.2*(1+relative_amplitude_variation*mrh_deviation)
-
- return Aw, min_rel_heat
-
- # *************************************************************************
- # *************************************************************************
-
- def get_single_ode_model_data(relative_amplitude_variation: float = 0.0):
-
- # define how the coefficients change
-
- Ria_deviation = random.random() - 0.5
-
- # define the (A, B, C and D) matrices
- # A: n*n
- # B: n*m
- # C: r*n
- # D: r*m
-
- Ci = 1.360*3600000
- Ria = (
- (1+relative_amplitude_variation*Ria_deviation)*5.31/3600000
- )
- Aw = 6.22
-
- min_rel_heat = 0.2
-
- x0 = np.array([20])
-
- return Ci, Ria, Aw, min_rel_heat, x0
-
- # *************************************************************************
- # *************************************************************************
-
- def get_multi_ode_model_data(
- self,
- relative_amplitude_variation: float = 0.0):
-
- # define how the coefficients change
-
- Rih_deviation = random.random() - 0.5
-
- Ria_deviation = random.random() - 0.5
-
- # define the (A, B, C and D) matrices
- # A: n*n
- # B: n*m
- # C: r*n
- # D: r*m
-
- # from Bacher and Madsen (2011): model TiTh
-
- Ci = 1.360*3600000 # original: 1.36 kWh/ºC
- Ch = 0.309*3600000 # original: 0.309 kWh/ºC
- Ria = (
- (1+relative_amplitude_variation*Ria_deviation)*5.31/3600000
- ) # original: 5.31 ºC/kWh
- Rih = (
- (1+relative_amplitude_variation*Rih_deviation)*0.639/3600000
- ) # original: 0.639 ºC/kWh
- Aw = 6.22 # original: 6.22 m2
-
- Pw = 5000 # 5 kW
-
- min_rel_heat = 0.2
-
- x0 = np.array([20, 20])
-
- return Ci, Ch, Ria, Rih, Aw, min_rel_heat, Pw, x0
-
- # *************************************************************************
- # *************************************************************************
-
- def stateless_model(self, Aw, min_rel_heat):
-
- # inputs: Ta, phi_s, phi_h above the minimum, phi_h status
- # outputs: solar irradiance, heat
-
- d = np.array([[0, Aw, 0, 0],
- [0, 0, (1-min_rel_heat), min_rel_heat]])
-
- return None, None, None, d
-
- # *************************************************************************
- # *************************************************************************
-
- def single_node_model(self, Ci, Ria, Aw, min_rel_heat, Pw):
-
- # states: Ti and Th
- # inputs: Ta, phi_s, phi_h above the minimum, phi_h status
- # outputs: solar irradiance, heat
-
- a = np.array([[-1/(Ria*Ci)]])
- b = np.array([[1/(Ci*Ria), Aw/Ci, Pw*(1-min_rel_heat)/Ci, Pw*min_rel_heat/Ci]])
- c = np.array([[0],[0]])
- d = np.array([[0, Aw, 0, 0],
- [0, 0, Pw*(1-min_rel_heat), Pw*min_rel_heat]])
-
- return a, b, c, d
-
- # *************************************************************************
- # *************************************************************************
-
- def two_node_model(self, Ci, Ch, Ria, Rih, Aw, min_rel_heat, Pw):
-
- # states: Ti and Th
- # inputs: Ta, phi_s, phi_h above the minimum, phi_h status
- # outputs: solar irradiance, heat
-
- a = np.array([[-(1/Rih+1/Ria)/Ci, 1/(Ci*Rih)],
- [1/(Ch*Rih), -1/(Ch*Rih)]])
- b = np.array([[1/(Ci*Ria), Aw/Ci, 0, 0],
- [0, 0, Pw*(1-min_rel_heat)/Ch, Pw*min_rel_heat/Ch]])
- c = np.array([[0, 0],
- [0, 0]])
- d = np.array([[0, Aw, 0, 0],
- [0, 0, Pw*(1-min_rel_heat), Pw*min_rel_heat]])
-
- return a, b, c, d
-
- # *************************************************************************
- # *************************************************************************
-
- def get_two_node_model_signals(self, number_samples):
-
- # signals
-
- # inputs:
- # 1) ambient temperature (real, can be fixed later)
- # 2) solar irradiation (real, can be fixed later)
- # 3) relative heat above minimum (nnr)
- # 4) heater status (binary)
-
- list_inputs = [
- sgn.FreeUnboundedSignal(number_samples),
- sgn.FreeUnboundedSignal(number_samples),
- sgn.NonNegativeRealSignal(number_samples),
- sgn.BinarySignal(number_samples)
- ]
-
- # states
- # 1) indoor temperature (real)
- # 2) heater temperature (real)
-
- list_states = [
- sgn.FreeUnboundedSignal(number_samples),
- sgn.FreeUnboundedSignal(number_samples)
- ]
-
- # outputs:
- # 1) solar gain (nnr)
- # 2) heat input (nnr)
-
- list_outputs = [
- sgn.NonNegativeRealSignal(number_samples),
- sgn.NonNegativeRealSignal(number_samples)
- ]
-
- return list_inputs, list_states, list_outputs
-
- # *************************************************************************
- # *************************************************************************
-
- def example_full_converter(self, time_step_durations: list):
-
- # number of samples
-
- number_time_steps = len(time_step_durations)
-
- # get the coefficients
-
- (Ci,
- Ch,
- Ria,
- Rih,
- Aw,
- min_rel_heat,
- Pw,
- x0) = self.get_multi_ode_model_data()
-
- # get the model
-
- (a,
- b,
- c,
- d) = self.two_node_model(Ci, Ch, Ria, Rih, Aw, min_rel_heat, Pw)
-
- # get the signals
-
- inputs, states, outputs = self.get_two_node_model_signals(
- number_time_steps)
-
- # create a dynamic system
-
- ds = dynsys.DynamicSystem(
- time_interval_durations=time_step_durations,
- A=a,
- B=b,
- C=c,
- D=d)
-
- # create a converter
-
- cvn1 = cvn.Converter(
- 'cvn1',
- sys=ds,
- initial_states=x0,
- turn_key_cost=3,
- inputs=inputs,
- states=states,
- outputs=outputs)
-
- # get the dictionaries
-
- (a_innk,
- b_inmk,
- c_irnk,
- d_irmk,
- e_x_ink,
- e_y_irk) = cvn1.matrix_dictionaries()
-
- # check the dicts
-
- #**************************************************************************
- #**************************************************************************
-
-#******************************************************************************
-#******************************************************************************
\ No newline at end of file
diff --git a/tests/test_data_finance.py b/tests/test_data_finance.py
index 521f009..e1b3af3 100644
--- a/tests/test_data_finance.py
+++ b/tests/test_data_finance.py
@@ -1,22 +1,219 @@
# imports
# local, internal
-
from src.topupopt.data.finance.invest import Investment, discount_factor, npv
-
from src.topupopt.data.finance.invest import salvage_value_linear_depreciation
-
from src.topupopt.data.finance.invest import present_salvage_value_annuity
-
from src.topupopt.data.finance.invest import salvage_value_annuity
+from src.topupopt.data.finance.utils import ArcInvestments
# local, external
-
import math
#******************************************************************************
#******************************************************************************
+class TestArcInvestments:
+
+ def test_object_creation(self):
+
+ discount_rate = 0.035
+ analysis_period = 20
+ discount_rates = tuple(
+ discount_rate for i in range(analysis_period)
+ )
+ cash_flows_npv = [
+ 3.673079208612225,
+ 1.5610827143687658,
+ 17,
+ 117,
+ 1842
+ ]
+ abs_cash_flows_npv = [
+ 1e-1,
+ 1e-1,
+ 1e-1,
+ 0.5,
+ 1
+ ]
+ investments = []
+
+ # limited longevity operation, does not go beyond planning horizon
+ cash_flow = 1
+ cash_flow_start = 1
+ myinv = Investment(discount_rates=discount_rates)
+ myinv.add_operational_cash_flows(
+ cash_flow=cash_flow,
+ start_period=cash_flow_start,
+ longevity=4
+ )
+ investments.append(myinv)
+
+ # limited longevity operation, goes beyond planning horizon
+ cash_flow = 1
+ cash_flow_start = 18
+ myinv = Investment(discount_rates=discount_rates)
+ myinv.add_operational_cash_flows(
+ cash_flow=cash_flow,
+ start_period=cash_flow_start,
+ longevity=4
+ )
+ investments.append(myinv)
+
+ # D&V omkostninger kundeanlæg: Sengeløse Skole
+ cash_flow = 1.2
+ cash_flow_start = 1
+ myinv = Investment(discount_rates=discount_rates)
+ myinv.add_operational_cash_flows(
+ cash_flow=cash_flow,
+ start_period=cash_flow_start,
+ longevity=None
+ )
+ investments.append(myinv)
+
+ # D&V omkostninger kundeanlæg: Større forbrugere
+ cash_flow = 6/10
+ myinv = Investment(discount_rates=discount_rates)
+ myinv.add_operational_cash_flows(
+ cash_flow=10*cash_flow,
+ start_period=1,
+ longevity=None
+ )
+ myinv.add_operational_cash_flows(
+ cash_flow=3*cash_flow,
+ start_period=2,
+ longevity=None
+ )
+ myinv.add_operational_cash_flows(
+ cash_flow=1*cash_flow,
+ start_period=3,
+ longevity=None
+ )
+ investments.append(myinv)
+
+ # D&V omkostninger kundeanlæg: Mindre forbrugere
+ cash_flow = 0.320
+ myinv = Investment(discount_rates=discount_rates)
+ myinv.add_operational_cash_flows(
+ cash_flow=197*cash_flow,
+ start_period=1,
+ longevity=None
+ )
+ myinv.add_operational_cash_flows(
+ cash_flow=(233-197)*cash_flow,
+ start_period=2,
+ longevity=None
+ )
+ myinv.add_operational_cash_flows(
+ cash_flow=(269-233)*cash_flow,
+ start_period=3,
+ longevity=None
+ )
+ myinv.add_operational_cash_flows(
+ cash_flow=(305-269)*cash_flow,
+ start_period=4,
+ longevity=None
+ )
+ myinv.add_operational_cash_flows(
+ cash_flow=(341-305)*cash_flow,
+ start_period=5,
+ longevity=None
+ )
+ myinv.add_operational_cash_flows(
+ cash_flow=(377-341)*cash_flow,
+ start_period=6,
+ longevity=None
+ )
+ myinv.add_operational_cash_flows(
+ cash_flow=(413-377)*cash_flow,
+ start_period=7,
+ longevity=None
+ )
+ myinv.add_operational_cash_flows(
+ cash_flow=(449-413)*cash_flow,
+ start_period=8,
+ longevity=None
+ )
+ myinv.add_operational_cash_flows(
+ cash_flow=(488-449)*cash_flow,
+ start_period=9,
+ longevity=None
+ )
+ investments.append(myinv)
+ # check
+ for inv, true_npv, _abs in zip(
+ investments,
+ cash_flows_npv,
+ abs_cash_flows_npv
+ ):
+ assert math.isclose(inv.net_present_value(), true_npv, abs_tol=_abs)
+ # create object
+ number_options = 5
+ static_loss = {
+ (h, q, k): 1+q*0.25+h*0.15+k*0.05
+ for h in range(number_options)
+ for q in range(2)
+ for k in range(3)
+ }
+ arc_invs = ArcInvestments(
+ investments=investments,
+ name='any',
+ efficiency=None,
+ efficiency_reverse=None,
+ capacity=tuple(1+o for o in range(number_options)),
+ #minimum_cost=tuple(1+o for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=False,
+ static_loss=static_loss,
+ validate=True
+ )
+ # check the costs
+ for _npv, true_npv, _abs in zip(
+ arc_invs.minimum_cost,
+ cash_flows_npv,
+ abs_cash_flows_npv
+ ):
+ assert math.isclose(_npv, true_npv, abs_tol=_abs)
+ # change something in the investments
+ for inv in arc_invs.investments:
+ inv.add_investment(
+ investment=1,
+ investment_period=0,
+ investment_longevity=10
+ )
+ # update the minimum costs
+ arc_invs.update_minimum_cost()
+ # make sure the comparison fails
+ for _npv, true_npv, _abs in zip(
+ arc_invs.minimum_cost,
+ cash_flows_npv,
+ abs_cash_flows_npv
+ ):
+ error_raised = False
+ try:
+ assert math.isclose(_npv, true_npv, abs_tol=_abs)
+ except AssertionError:
+ error_raised = True
+ assert error_raised
+ # new true values
+ new_true_npv = [
+ 4.673079208612225,
+ 2.561082714368766,
+ 18.054883962342764,
+ 117.50524073646935,
+ 1843.1804103901486
+ ]
+ # print([inv.net_present_value() for inv in arc_invs])
+ for _npv, true_npv, _abs in zip(
+ arc_invs.minimum_cost,
+ new_true_npv,
+ abs_cash_flows_npv
+ ):
+ assert math.isclose(_npv, true_npv, abs_tol=_abs)
+
+# *****************************************************************************
+# *****************************************************************************
+
class TestDataFinance:
# TODO: make sure that all methods work with variable discount rates
diff --git a/tests/test_dhn.py b/tests/test_dhn.py
index 893128a..96e4263 100644
--- a/tests/test_dhn.py
+++ b/tests/test_dhn.py
@@ -11,21 +11,17 @@ from math import inf, isclose
# local, internal
# topupopt
-
-# import src.topupopt.problems.esipp as tuo
-
from src.topupopt.data.dhn.network import PipeTrenchOptions, ExistingPipeTrench
-# from src.topupopt.data.dhn.network import
+from src.topupopt.data.dhn.network import PipeTrenchInvestments
+from src.topupopt.data.finance.invest import Investment
# topupheat
-
from topupheat.pipes.single import StandardisedPipe, StandardisedPipeDatabase
-
import topupheat.pipes.trenches as trenches
-from topupheat.common.fluids import FluidDatabase#, Fluid
+from topupheat.common.fluids import FluidDatabase
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
class TestDistrictHeatingNetwork:
@@ -147,6 +143,9 @@ class TestDistrictHeatingNetwork:
# *********************************************************************
+ # *************************************************************************
+ # *************************************************************************
+
def test_creating_single_arcs(self):
# fluid data
@@ -283,8 +282,156 @@ class TestDistrictHeatingNetwork:
outdoor_temperature for i in range(number_steps)
]
)
+
+ # *************************************************************************
+ # *************************************************************************
+ def test_creating_single_arcs_investment(self):
+
+ # fluid data
+ waterdata_file = 'tests/data/incropera2006_saturated_water.csv'
+ phase = FluidDatabase.fluid_LIQUID
+ fluid_db = FluidDatabase(
+ fluid='fluid',
+ phase=phase,
+ source=waterdata_file
+ )
+
+ singlepipedata_files = ['tests/data/isoplus_singlepipes_s1.csv']
+ pipedb = StandardisedPipeDatabase(source=singlepipedata_files)
+ pipe = StandardisedPipe(
+ pipe_tuple=pipedb.pipe_tuples[0],
+ #e_eff=pipe_e_eff,
+ #sp=pipe_specific_price,
+ db=pipedb)
+
+ # network details
+ supply_temperature = 85+273.15
+ return_temperature = 45+273.15
+ pressure = 1e5
+ # trench
+ pipe_distance = 0.52 # m
+ pipe_depth = 0.66 # m
+ # environmental
+ outdoor_temperature = 6+273.15 # K
+ h_gs = inf # 14.6 # W/m2K
+ soil_k = 1.5 # W/mK
+ # more information
+ max_specific_pressure_loss = 100 # Pa/m
+
+ mytrench = trenches.SupplyReturnPipeTrench(
+ pipe_center_depth=pipe_depth,
+ pipe_center_distance=pipe_distance,
+ fluid_db=fluid_db,
+ phase=phase,
+ pressure=pressure,
+ supply_temperature=supply_temperature,
+ return_temperature=return_temperature,
+ max_specific_pressure_loss=max_specific_pressure_loss,
+ supply_pipe=pipe
+ )
+
+ # investments
+ number_periods = 20
+ discount_rate = 0.035
+ discount_rates = tuple([discount_rate for p in range(number_periods)])
+ inv = Investment(discount_rates=discount_rates)
+
+ # PipeTrenchOptions
+ myarcs = PipeTrenchInvestments(
+ trench=mytrench,
+ name='hellotrench',
+ length=50,
+ investments=(inv,),
+ )
+
+ # add static loss scenario
+ myarcs.set_static_losses(
+ scenario_key='scenario0',
+ ground_thermal_conductivity=soil_k,
+ ground_air_heat_transfer_coefficient=h_gs,
+ time_interval_duration=3600,
+ temperature_surroundings=outdoor_temperature
+ )
+ # add another static loss scenario
+ myarcs.set_static_losses(
+ scenario_key='scenario1',
+ ground_thermal_conductivity=soil_k+1,
+ ground_air_heat_transfer_coefficient=h_gs+1,
+ time_interval_duration=3600+100,
+ temperature_surroundings=outdoor_temperature+1
+ )
+ number_steps = 3
+ myarcs.set_static_losses(
+ scenario_key='scenario2',
+ ground_thermal_conductivity=[soil_k for i in range(number_steps)
+ ],
+ ground_air_heat_transfer_coefficient=[
+ h_gs for i in range(number_steps)
+ ],
+ time_interval_duration=[3600 for i in range(number_steps)],
+ temperature_surroundings=[
+ outdoor_temperature for i in range(number_steps)
+ ]
+ )
+
+ # test arcs
+
+ n = myarcs.number_options()
+ assert not myarcs.has_been_selected()
+ assert len(myarcs.capacity) == n
+ assert len(myarcs.minimum_cost) == n
+ assert isinstance(myarcs.specific_capacity_cost, Real)
+ assert type(myarcs.efficiency) == type(None)
+ assert type(myarcs.efficiency_reverse) == type(None)
+ assert type(myarcs.static_loss) == dict
+ assert len(myarcs.static_loss) != 0
+ for (h, q, k), sl in myarcs.static_loss.items():
+ assert isinstance(sl, Real)
+ assert sl >= 0
+
+ # redefine the capacity
+ capacity = tuple(myarcs.capacity)
+ myarcs.set_capacity(
+ max_specific_pressure_loss=max_specific_pressure_loss+100
+ )
+ assert len(capacity) == len(myarcs.capacity)
+ for _c1, _c2 in zip(capacity, myarcs.capacity):
+ assert _c1 != _c2
+
+ # redefine the minimum costs
+ min_cost = tuple(myarcs.minimum_cost)
+ myarcs.set_minimum_cost(minimum_cost=[_mc+1 for _mc in min_cost])
+ assert len(min_cost) == len(myarcs.minimum_cost)
+ for _mc1, _mc2 in zip(min_cost, myarcs.minimum_cost):
+ assert _mc1+1 == _mc2
+
# *********************************************************************
+
+ # create arcs object with multiple static loss values as a first case
+
+ # PipeTrenchOptions
+ myarcs = PipeTrenchOptions(
+ trench=mytrench,
+ name='hellotrench',
+ length=50
+ )
+
+ number_steps = 3
+ myarcs.set_static_losses(
+ scenario_key='scenario2',
+ ground_thermal_conductivity=[soil_k for i in range(number_steps)],
+ ground_air_heat_transfer_coefficient=[
+ h_gs for i in range(number_steps)
+ ],
+ time_interval_duration=[3600 for i in range(number_steps)],
+ temperature_surroundings=[
+ outdoor_temperature for i in range(number_steps)
+ ]
+ )
+
+ # *************************************************************************
+ # *************************************************************************
def test_creating_multiple_arcs(self):
@@ -444,7 +591,8 @@ class TestDistrictHeatingNetwork:
]
)
- # *********************************************************************
+ # *************************************************************************
+ # *************************************************************************
# *****************************************************************************
# *****************************************************************************
diff --git a/tests/test_dhn_utils.py b/tests/test_dhn_utils.py
index b2e5139..147428b 100644
--- a/tests/test_dhn_utils.py
+++ b/tests/test_dhn_utils.py
@@ -18,6 +18,220 @@ from topupheat.common.fluids import FluidDatabase#, Fluid
# *****************************************************************************
class TestDistrictHeatingNetworkUtils:
+
+ # *************************************************************************
+ # *************************************************************************
+
+ def test_cost_pipes_single_arc(self):
+
+ # fluid data
+ waterdata_file = 'tests/data/incropera2006_saturated_water.csv'
+ phase = FluidDatabase.fluid_LIQUID
+ fluid_db = FluidDatabase(
+ fluid='fluid',
+ phase=phase,
+ source=waterdata_file
+ )
+
+ singlepipedata_files = ['tests/data/isoplus_singlepipes_s1.csv']
+ pipedb = StandardisedPipeDatabase(source=singlepipedata_files)
+ pipe = StandardisedPipe(
+ pipe_tuple=pipedb.pipe_tuples[0],
+ #e_eff=pipe_e_eff,
+ #sp=pipe_specific_price,
+ db=pipedb)
+
+ # network details
+ supply_temperature = 85+273.15
+ return_temperature = 45+273.15
+ pressure = 1e5
+ # trench
+ pipe_distance = 0.52 # m
+ pipe_depth = 0.66 # m
+ # environmental
+ outdoor_temperature = 6+273.15 # K
+ h_gs = inf # 14.6 # W/m2K
+ soil_k = 1.5 # W/mK
+ # more information
+ max_specific_pressure_loss = 100 # Pa/m
+
+ mytrench = trenches.SupplyReturnPipeTrench(
+ pipe_center_depth=pipe_depth,
+ pipe_center_distance=pipe_distance,
+ fluid_db=fluid_db,
+ phase=phase,
+ pressure=pressure,
+ supply_temperature=supply_temperature,
+ return_temperature=return_temperature,
+ max_specific_pressure_loss=max_specific_pressure_loss,
+ supply_pipe=pipe
+ )
+
+ # create arcs object with multiple static loss values as a first case
+
+ trench_length = 50
+
+ # PipeTrenchOptions
+ myarcs = PipeTrenchOptions(
+ trench=mytrench,
+ name='hellotrench',
+ length=trench_length
+ )
+
+ number_steps = 3
+ myarcs.set_static_losses(
+ scenario_key='scenario2',
+ ground_thermal_conductivity=[soil_k for i in range(number_steps)],
+ ground_air_heat_transfer_coefficient=[
+ h_gs for i in range(number_steps)
+ ],
+ time_interval_duration=[3600 for i in range(number_steps)],
+ temperature_surroundings=[
+ outdoor_temperature for i in range(number_steps)
+ ]
+ )
+
+ mypipecosts = utils.cost_pipes(mytrench, trench_length)
+ assert mypipecosts == (50.0,)
+
+ # unrecognised input: using a string for the trench length
+ error_raised = False
+ try:
+ mypipecosts = utils.cost_pipes(mytrench, '50')
+ except ValueError:
+ error_raised = True
+ assert error_raised
+
+ # *************************************************************************
+ # *************************************************************************
+
+ def test_cost_pipes_multiple_arcs(self):
+
+ # fluid data
+ waterdata_file = 'tests/data/incropera2006_saturated_water.csv'
+ phase = FluidDatabase.fluid_LIQUID
+ fluid_db = FluidDatabase(
+ fluid='fluid',
+ phase=phase,
+ source=waterdata_file
+ )
+
+ singlepipedata_files = ['tests/data/isoplus_singlepipes_s1.csv']
+ pipedb = StandardisedPipeDatabase(source=singlepipedata_files)
+ pipe = StandardisedPipe(
+ pipe_tuple=pipedb.pipe_tuples[0],
+ db=pipedb)
+
+ # network details
+ supply_temperature = 85+273.15
+ return_temperature = 45+273.15
+ pressure = 1e5
+ # trench
+ pipe_distance = 0.52 # m
+ pipe_depth = 0.66 # m
+ # environmental
+ outdoor_temperature = 6+273.15 # K
+ h_gs = inf # 14.6 # W/m2K
+ soil_k = 1.5 # W/mK
+ # more information
+ max_specific_pressure_loss = 100 # Pa/m
+ number_options = 2
+
+ mytrench = trenches.SupplyReturnPipeTrench(
+ pipe_center_depth=[pipe_depth for i in range(number_options)],
+ pipe_center_distance=[
+ pipe_distance for i in range(number_options)
+ ],
+ fluid_db=fluid_db,
+ phase=phase,
+ pressure=[pressure for i in range(number_options)],
+ supply_temperature=[
+ supply_temperature for i in range(number_options)
+ ],
+ return_temperature=[
+ return_temperature for i in range(number_options)
+ ],
+ max_specific_pressure_loss=[
+ max_specific_pressure_loss for i in range(number_options)
+ ],
+ supply_pipe=[pipe for i in range(number_options)]
+ )
+
+ # PipeTrenchOptions
+ myarcs = PipeTrenchOptions(
+ trench=mytrench,
+ name='hellotrench',
+ length=50
+ )
+
+ # add static loss scenario
+ myarcs.set_static_losses(
+ scenario_key='scenario0',
+ ground_thermal_conductivity=soil_k,
+ ground_air_heat_transfer_coefficient=h_gs,
+ time_interval_duration=3600,
+ temperature_surroundings=outdoor_temperature
+ )
+ # add another static loss scenario
+ myarcs.set_static_losses(
+ scenario_key='scenario1',
+ ground_thermal_conductivity=soil_k+1,
+ ground_air_heat_transfer_coefficient=h_gs+1,
+ time_interval_duration=3600+100,
+ temperature_surroundings=outdoor_temperature+1
+ )
+ # add static loss scenario
+ number_steps = 3
+ myarcs.set_static_losses(
+ scenario_key='scenario2',
+ ground_thermal_conductivity=[soil_k for i in range(number_steps)
+ ],
+ ground_air_heat_transfer_coefficient=[
+ h_gs for i in range(number_steps)
+ ],
+ time_interval_duration=[3600 for i in range(number_steps)],
+ temperature_surroundings=[
+ outdoor_temperature for i in range(number_steps)
+ ]
+ )
+ trench_length = 50
+
+ # PipeTrenchOptions
+ myarcs = PipeTrenchOptions(
+ trench=mytrench,
+ name='hellotrench',
+ length=trench_length
+ )
+
+ number_steps = 3
+ myarcs.set_static_losses(
+ scenario_key='scenario2',
+ ground_thermal_conductivity=[soil_k for i in range(number_steps)
+ ],
+ ground_air_heat_transfer_coefficient=[
+ h_gs for i in range(number_steps)
+ ],
+ time_interval_duration=[3600 for i in range(number_steps)],
+ temperature_surroundings=[
+ outdoor_temperature for i in range(number_steps)
+ ]
+ )
+
+ mypipecosts = utils.cost_pipes(mytrench, trench_length)
+ assert mypipecosts == (100.0, 100.0)
+ mypipecosts = utils.cost_pipes(mytrench, (trench_length, trench_length))
+ assert mypipecosts == (100.0, 100.0)
+
+ # unrecognised input: using a list for the trench lengths
+ error_raised = False
+ try:
+ mypipecosts = utils.cost_pipes(mytrench, [trench_length])
+ except ValueError:
+ error_raised = True
+ assert error_raised
+
+ # *************************************************************************
+ # *************************************************************************
def test_plotting_heating_demand(self):
@@ -272,10 +486,12 @@ class TestDistrictHeatingNetworkUtils:
# *********************************************************************
- out = utils.summarise_network_by_arc_technology(mynet, False)
+ out = utils.summarise_network_by_pipe_technology(mynet, False)
assert 'DN20' in out and out['DN20'] == 25
assert 'DN32' in out and out['DN32'] == 50
+
+ utils.summarise_network_by_pipe_technology(mynet, True)
# *************************************************************************
# *************************************************************************
@@ -427,13 +643,18 @@ class TestDistrictHeatingNetworkUtils:
# *********************************************************************
- utils.summarise_network_by_arc_technology(network, False)
+ utils.summarise_network_by_pipe_technology(network, False)
utils.plot_network_layout(
network=network,
include_basemap=False)
+
+ utils.plot_network_layout(
+ network=network,
+ include_basemap=True)
# *************************************************************************
- # *************************************************************************
+ # *************************************************************************
+
#******************************************************************************
#******************************************************************************
\ No newline at end of file
diff --git a/tests/examples_converter.py b/tests/test_esipp_converter.py
similarity index 88%
rename from tests/examples_converter.py
rename to tests/test_esipp_converter.py
index dfb7b08..2e1f9ac 100644
--- a/tests/examples_converter.py
+++ b/tests/test_esipp_converter.py
@@ -19,7 +19,7 @@ import src.topupopt.problems.esipp.signal as sgn
#******************************************************************************
#******************************************************************************
-def examples():
+class TestConverter():
#**************************************************************************
#**************************************************************************
@@ -47,10 +47,23 @@ def examples():
# test creating a converter based on a multiple ODE system with 1 output
# test creating a converter based on a multiple ODE multi-output system
- time_step_durations = [1,1,1,1]
- example_full_converter(time_step_durations)
- time_step_durations = [1,1.5,0.5,1]
- example_full_converter(time_step_durations)
+
+
+ #**************************************************************************
+ #**************************************************************************
+
+ def test_full_converter_regular(self):
+
+ time_step_durations = [1,1,1,1]
+ method_full_converter(time_step_durations)
+
+ #**************************************************************************
+ #**************************************************************************
+
+ def test_full_converter_irregular(self):
+
+ time_step_durations = [1,1.5,0.5,1]
+ method_full_converter(time_step_durations)
#**************************************************************************
#**************************************************************************
@@ -95,8 +108,8 @@ def get_single_ode_model_data(relative_amplitude_variation: float = 0.0):
return Ci, Ria, Aw, min_rel_heat, x0
-#******************************************************************************
-#******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
def get_multi_ode_model_data(relative_amplitude_variation: float = 0.0):
@@ -225,26 +238,21 @@ def get_two_node_model_signals(number_samples):
#******************************************************************************
#******************************************************************************
-def example_full_converter(time_step_durations: list):
+def method_full_converter(time_step_durations: list):
# number of samples
-
number_time_steps = len(time_step_durations)
# get the coefficients
-
Ci, Ch, Ria, Rih, Aw, min_rel_heat, Pw, x0 = get_multi_ode_model_data()
# get the model
-
a, b, c, d = two_node_model(Ci, Ch, Ria, Rih, Aw, min_rel_heat, Pw)
# get the signals
-
inputs, states, outputs = get_two_node_model_signals(number_time_steps)
# create a dynamic system
-
ds = dynsys.DynamicSystem(
time_interval_durations=time_step_durations,
A=a,
@@ -253,7 +261,6 @@ def example_full_converter(time_step_durations: list):
D=d)
# create a converter
-
cvn1 = cvn.Converter(
'cvn1',
sys=ds,
@@ -264,11 +271,14 @@ def example_full_converter(time_step_durations: list):
outputs=outputs)
# get the dictionaries
-
- a_innk, b_inmk, c_irnk, d_irmk, e_x_ink, e_y_irk = cvn1.matrix_dictionaries()
-
- # check the dicts
-
+ (a_innk,
+ b_inmk,
+ c_irnk,
+ d_irmk,
+ e_x_ink,
+ e_y_irk) = cvn1.matrix_dictionaries()
+
+ # TODO: check the dicts
#******************************************************************************
#******************************************************************************
\ No newline at end of file
diff --git a/tests/test_esipp_dynsys.py b/tests/test_esipp_dynsys.py
new file mode 100644
index 0000000..756e1f7
--- /dev/null
+++ b/tests/test_esipp_dynsys.py
@@ -0,0 +1,2231 @@
+# imports
+
+# standard
+import random
+
+# local, external
+import numpy as np
+from scipy.signal import StateSpace, lsim, dlsim
+
+# local, internal
+import src.topupopt.problems.esipp.dynsys as dynsys
+
+# *****************************************************************************
+# *****************************************************************************
+
+class TestDynsys():
+
+ # *************************************************************************
+ # *************************************************************************
+
+# seed_number = random.randint(1,int(1e5))
+
+# print_outputs = True
+
+# # with states and outputs
+
+# # test multi-ODE, multi-output dynamic systems while integrating outputs
+
+# examples_dynsys_multiode_multiout(True, print_outputs, seed_number)
+
+# # test multi-ODE, multi-output dynamic systems without integrating outputs
+
+# examples_dynsys_multiode_multiout(False, print_outputs, seed_number)
+
+# # test single ODE, multi-output dynamic systems while integrating outputs
+
+# examples_dynsys_singleode_multiout(True, print_outputs, seed_number)
+
+# # test single ODE, multi-output dynamic systems without integrating outputs
+
+# examples_dynsys_singleode_multiout(False, print_outputs, seed_number)
+
+# # test multi-ODE, single-output dynamic systems while integrating outputs
+
+# examples_dynsys_multiode_multiout(True, print_outputs, seed_number, 1)
+
+# # test multi-ODE, single-output dynamic systems without integrating outputs
+
+# examples_dynsys_multiode_multiout(False, print_outputs, seed_number, 1)
+
+# # test single-ODE, single-output dynamic systems while integrating outputs
+
+# examples_dynsys_singleode_multiout(True, print_outputs, seed_number, 1)
+
+# # test single-ODE, single-output dynamic systems without integrating outputs
+
+# examples_dynsys_singleode_multiout(False, print_outputs, seed_number, 1)
+
+# # *************************************************************************
+
+# # outputless
+
+# # test single-ODE, outputless dynamic systems while integrating outputs
+
+# examples_dynsys_singleode_multiout(True, print_outputs, seed_number, 0)
+
+# # test multi-ODE, outputless dynamic systems while integrating outputs
+
+# examples_dynsys_multiode_multiout(True, print_outputs, seed_number, 0)
+
+# # test single-ODE, outputless dynamic systems without integrating outputs
+
+# examples_dynsys_singleode_multiout(False, print_outputs, seed_number, 0)
+
+# # test multi-ODE, outputless dynamic systems without integrating outputs
+
+# examples_dynsys_multiode_multiout(False, print_outputs, seed_number, 0)
+
+# # outputless system via dynsys subclass
+
+# example_outputless_system_object()
+
+# # *************************************************************************
+
+# # stateless
+
+# # test stateless, single-output dynamic systems while integrating outputs
+
+# examples_dynsys_stateless_multiout(True, print_outputs, seed_number, 1)
+
+# # test stateless, multi-output dynamic systems without integrating outputs
+
+# examples_dynsys_stateless_multiout(False, print_outputs, seed_number, 2)
+
+# # stateless system via dynsys subclass
+
+# example_stateless_system_object(True)
+# example_stateless_system_object(False)
+
+# # *************************************************************************
+# # *************************************************************************
+
+# # trigger errors
+
+# # test stateless, outputless dynamic systems while integrating outputs
+
+# number_errors = 0
+
+# try:
+# examples_dynsys_stateless_multiout(True, False, seed_number, 0)
+# except Exception:
+# number_errors += 1
+
+# assert number_errors == 1
+
+# # test stateless, outputless dynamic systems without integrating outputs
+
+# number_errors = 0
+
+# try:
+# examples_dynsys_stateless_multiout(False, False, seed_number, 0)
+# except Exception:
+# number_errors += 1
+
+# assert number_errors == 1
+
+# # test negative time duration
+
+# example_incorrect_time_step_durations()
+
+# # test unrecognised matrix formats
+
+# example_unrecognised_matrix_formats()
+
+# # different matrix sizes for the same problem, all other things being equal
+
+# example_varying_matrix_sizes(True)
+# example_varying_matrix_sizes(False)
+
+# # test multiple A matrices and multiple non-matching time intervals
+
+# example_nonmatching_time_steps_and_matrices()
+
+# # test non-square A matrices
+
+# example_nonsquare_A_matrices()
+
+# # test incompatible A and B matrices (different number of rows)
+
+# example_incompatible_AB_matrices()
+
+# # test incompatible C and D matrices (different number of rows)
+
+# example_incompatible_CD_matrices()
+
+# # test incompatible A and C matrices (different number of columns)
+
+# example_incompatible_AC_matrices()
+
+# # test incompatible B and D matrices (different number of columns)
+
+# example_incompatible_BD_matrices()
+
+# # trigger incorrect input signal format error when simulating
+
+# example_single_time_step_model_incorrect_inputs()
+
+# # TODO: test only some matrices as being time invariant
+
+ # *************************************************************************
+ # *************************************************************************
+
+ def test_incorrect_time_step_durations(self):
+
+ integrate_outputs = True
+
+ # *********************************************************************
+
+ # negative time step duration
+
+ # regular time steps (as an int)
+
+ time_step_durations = [-1]
+
+ # get the model
+
+ Aw, min_rel_heat = get_stateless_model_data()
+
+ _, _, _, D = stateless_model(Aw, min_rel_heat)
+
+ number_errors = 0
+ try:
+ _ = dynsys.StatelessSystem(
+ time_interval_durations=time_step_durations[0],
+ D=D,
+ integrate_outputs=integrate_outputs)
+ except ValueError:
+ number_errors += 1
+
+ assert number_errors == 1
+
+ # *********************************************************************
+
+ # no time step duration (empty list)
+
+ time_step_durations = []
+
+ # get the model
+
+ Aw, min_rel_heat = get_stateless_model_data()
+
+ _, _, _, D = stateless_model(Aw, min_rel_heat)
+
+ number_errors = 0
+ try:
+ _ = dynsys.StatelessSystem(
+ time_interval_durations=time_step_durations,
+ D=D,
+ integrate_outputs=integrate_outputs)
+ except ValueError:
+ number_errors += 1
+
+ assert number_errors == 1
+
+ # *********************************************************************
+
+ # time step duration list with negative or zero time step durations
+
+ time_step_durations = [-1, 3, 0]
+
+ # get the model
+
+ Aw, min_rel_heat = get_stateless_model_data()
+
+ _, _, _, D = stateless_model(Aw, min_rel_heat)
+
+ number_errors = 0
+ try:
+ _ = dynsys.StatelessSystem(
+ time_interval_durations=time_step_durations,
+ D=D,
+ integrate_outputs=integrate_outputs)
+ except ValueError:
+ number_errors += 1
+
+ assert number_errors == 1
+
+ # *********************************************************************
+
+ # time step duration list with non-numeric types
+
+ time_step_durations = [None, 3, 3]
+
+ # get the model
+
+ Aw, min_rel_heat = get_stateless_model_data()
+
+ _, _, _, D = stateless_model(Aw, min_rel_heat)
+
+ number_errors = 0
+ try:
+ _ = dynsys.StatelessSystem(
+ time_interval_durations=time_step_durations,
+ D=D,
+ integrate_outputs=integrate_outputs)
+ except TypeError:
+ number_errors += 1
+
+ assert number_errors == 1
+
+ # *************************************************************************
+ # *************************************************************************
+
+ def test_single_time_step_model_incorrect_inputs(self):
+
+ integrate_y = True
+
+ # regular time steps (as an int)
+
+ time_step_durations = [1]
+
+ # define a sequence of time steps
+
+ t = np.array([sum(time_step_durations[0:i])
+ for i in range(len(time_step_durations)+1)])
+
+ # define the inputs
+
+ list_inputs = [-5, 50, 0.1, 1, 5] # extra input signal to force error
+
+ U = np.array([[input_i for dt in t] # extra time step to force special case
+ for input_i in list_inputs])
+
+ # U = np.array([[input_i for _ in range(len(time_step_durations))]
+ # for input_i in list_inputs])
+
+ # get the model
+
+ Ci, Ch, Ria, Rih, Aw, min_rel_heat, x0 = get_multi_ode_model_data()
+
+ A, B, C, D = two_node_model(Ci, Ch, Ria, Rih, Aw, min_rel_heat)
+
+ ds = dynsys.DynamicSystem(
+ time_interval_durations=time_step_durations,
+ A=A,
+ B=B,
+ C=C,
+ D=D,
+ integrate_outputs=integrate_y)
+
+ # define the initial conditions
+
+ number_errors = 0
+
+ try:
+ X, Y = ds.simulate(U, x0)
+ except ValueError:
+ number_errors += 1
+
+ assert number_errors == 1
+
+ # *************************************************************************
+ # *************************************************************************
+
+ def test_outputless_system_object(self):
+
+ # regular time steps
+ time_step_durations = [1, 1, 1, 1]
+
+ # define a sequence of time steps
+ t = np.array([sum(time_step_durations[0:i])
+ for i in range(len(time_step_durations)+1)])
+
+ # define the inputs
+ list_inputs = [-5, 50, 0.1, 1]
+ U = np.array([[input_i for _ in range(len(time_step_durations))]
+ for input_i in list_inputs])
+
+ # get the model
+ Ci, Ch, Ria, Rih, Aw, min_rel_heat, x0 = get_multi_ode_model_data()
+ A, B, _, _ = two_node_model(Ci, Ch, Ria, Rih, Aw, min_rel_heat)
+ ds = dynsys.OutputlessSystem(
+ time_interval_durations=time_step_durations,
+ A=A,
+ B=B)
+
+ # define the initial conditions
+ X, Y = ds.simulate(U, x0)
+ assert Y == None
+ assert isinstance(X, np.ndarray)
+
+ # *********************************************************************
+
+ # regular time steps (as an int)
+ time_step_durations = [1]
+
+ # define a sequence of time steps
+ t = np.array([sum(time_step_durations[0:i])
+ for i in range(len(time_step_durations)+1)])
+
+ # define the inputs
+ list_inputs = [-5, 50, 0.1, 1]
+ U = np.array([[input_i for _ in range(len(time_step_durations))]
+ for input_i in list_inputs])
+
+ # get the model
+ Ci, Ch, Ria, Rih, Aw, min_rel_heat, x0 = get_multi_ode_model_data()
+ A, B, _, _ = two_node_model(Ci, Ch, Ria, Rih, Aw, min_rel_heat)
+
+ ds = dynsys.OutputlessSystem(
+ time_interval_durations=time_step_durations,
+ A=A,
+ B=B)
+
+ # define the initial conditions
+ X, Y = ds.simulate(U, x0)
+ assert Y == None
+ assert isinstance(X, np.ndarray)
+
+ # *********************************************************************
+
+ # irregular time steps
+ time_step_durations = [1, 1.5, 0.5, 1]
+
+ # define a sequence of time steps
+ t = np.array([sum(time_step_durations[0:i])
+ for i in range(len(time_step_durations)+1)])
+
+ # define the inputs
+ list_inputs = [-5, 50, 0.1, 1]
+ U = np.array([[input_i for _ in range(len(time_step_durations))]
+ for input_i in list_inputs])
+
+ # get the model
+ Ci, Ch, Ria, Rih, Aw, min_rel_heat, x0 = get_multi_ode_model_data()
+ A, B, _, _ = two_node_model(Ci, Ch, Ria, Rih, Aw, min_rel_heat)
+ ds = dynsys.OutputlessSystem(
+ time_interval_durations=time_step_durations,
+ A=A,
+ B=B)
+
+ # define the initial conditions
+ X, Y = ds.simulate(U, x0)
+ assert Y == None
+ assert isinstance(X, np.ndarray)
+
+ # *************************************************************************
+ # *************************************************************************
+
+ def test_stateless_system_object_integration(self):
+
+ integrate_outputs = True
+ method_stateless_system_object(integrate_outputs)
+
+ # *************************************************************************
+ # *************************************************************************
+
+ def test_stateless_system_object_no_integration(self):
+
+ integrate_outputs = False
+ method_stateless_system_object(integrate_outputs)
+
+ # *************************************************************************
+ # *************************************************************************
+
+ def test_incompatible_BD_matrices(self):
+
+ integrate_y = True
+
+ # B and D must have the same number of columns
+
+ time_step_durations = [1, 1, 1]
+
+ A = np.array([[4, 9], [1, 3]])
+
+ B = np.array([[5, 6, 1], [7, 2, 1]])
+
+ C = np.array([[4, 6]])
+
+ D = np.array([[3, 6]])
+
+ # test
+
+ number_errors = 0
+
+ try:
+ _ = dynsys.DynamicSystem(
+ time_interval_durations=time_step_durations,
+ A=A,
+ B=B,
+ C=C,
+ D=D,
+ integrate_outputs=integrate_y)
+ except ValueError:
+ number_errors += 1
+
+ assert number_errors == 1
+
+ # *************************************************************************
+ # *************************************************************************
+
+ def test_incompatible_AC_matrices(self):
+
+ integrate_y = True
+
+ # A and C must have the same number of columns
+
+ time_step_durations = [1, 1, 1]
+
+ A = np.array([[4, 9], [1, 3]])
+
+ B = np.array([[5, 6, 1], [7, 2, 1]])
+
+ C = np.array([[4]])
+
+ D = np.array([[3, 6, 3]])
+
+ # test
+
+ number_errors = 0
+
+ try:
+ _ = dynsys.DynamicSystem(
+ time_interval_durations=time_step_durations,
+ A=A,
+ B=B,
+ C=C,
+ D=D,
+ integrate_outputs=integrate_y)
+ except ValueError:
+ number_errors += 1
+
+ assert number_errors == 1
+
+ # *****************************************************************************
+ # *****************************************************************************
+
+ def test_incompatible_CD_matrices(self):
+
+ integrate_y = True
+
+ # C and D must have the same number of rows
+
+ time_step_durations = [1, 1, 8]
+
+ A = np.array([[4, 9], [1, 3]])
+
+ B = np.array([[5, 6, 1], [7, 2, 1]])
+
+ C = np.array([[4, 9], [1, 3]])
+
+ D = np.array([[3, 6, 3]])
+
+ # test
+
+ number_errors = 0
+
+ try:
+ _ = dynsys.DynamicSystem(
+ time_interval_durations=time_step_durations,
+ A=A,
+ B=B,
+ C=C,
+ D=D,
+ integrate_outputs=integrate_y)
+ except ValueError:
+ number_errors += 1
+
+ assert number_errors == 1
+
+ # *************************************************************************
+ # *************************************************************************
+
+ def test_incompatible_AB_matrices(self):
+
+ integrate_y = True
+
+ # A and B must have the same number of rows
+
+ time_step_durations = [1, 1, 8]
+
+ A = np.array([[4, 9],[1, 3]])
+
+ B = np.array([[3, 6, 3]])
+
+ C = np.array([[2, 4]])
+
+ D = np.array([[5, 6, 1]])
+
+ # test A
+
+ number_errors = 0
+
+ try:
+ _ = dynsys.DynamicSystem(
+ time_interval_durations=time_step_durations,
+ A=A,
+ B=B,
+ C=C,
+ D=D,
+ integrate_outputs=integrate_y)
+ except ValueError:
+ number_errors += 1
+
+ assert number_errors == 1
+
+ # *************************************************************************
+ # *************************************************************************
+
+ def test_nonsquare_A_matrices(self):
+
+ integrate_y = True
+
+ # Single non-square A matrix
+
+ time_step_durations = [1, 1, 8]
+
+ A = np.array([[4, 9]])
+
+ B = np.array([[3]])
+
+ C = np.array([[2]])
+
+ D = np.array([[5]])
+
+ # test A
+
+ number_errors = 0
+
+ try:
+ _ = dynsys.DynamicSystem(
+ time_interval_durations=time_step_durations,
+ A=A,
+ B=B,
+ C=C,
+ D=D,
+ integrate_outputs=integrate_y)
+ except ValueError:
+ number_errors += 1
+
+ assert number_errors == 1
+
+ # # multiple matrices, at least one of which is non-square
+
+ # A = [np.array([[4]]), np.array([[1, 9]])]
+
+ # B = [np.array([[3]]), np.array([[4]])]
+
+ # C = [np.array([[2]]), np.array([[2]])]
+
+ # D = [np.array([[8]]), np.array([[7]])]
+
+ # # test A
+
+ # number_errors = 0
+
+ # try:
+ # _ = dynsys.DynamicSystem(
+ # time_interval_durations=time_step_durations,
+ # A=A,
+ # B=B,
+ # C=C,
+ # D=D,
+ # integrate_outputs=integrate_y)
+ # except ValueError:
+ # number_errors += 1
+
+ # assert number_errors == 1
+
+ # *************************************************************************
+ # *************************************************************************
+
+ def test_nonmatching_time_steps_and_matrices(self):
+
+ integrate_y = True
+
+ # multiple matrices and time intervals, but more time steps than matrices
+
+ time_step_durations = [1, 1, 8]
+
+ A = [np.array([[4]]),np.array([[1]])]
+
+ B = [np.array([[3]]),np.array([[4]]),np.array([[4]])]
+
+ C = [np.array([[2]]),np.array([[2]]),np.array([[2]])]
+
+ D = [np.array([[8]]),np.array([[7]]),np.array([[7]])]
+
+ # test A
+
+ number_errors = 0
+
+ try:
+ _ = dynsys.DynamicSystem(
+ time_interval_durations=time_step_durations,
+ A=A,
+ B=B,
+ C=C,
+ D=D,
+ integrate_outputs=integrate_y)
+ except ValueError:
+ number_errors += 1
+
+ assert number_errors == 1
+
+ # test B
+
+ number_errors = 0
+
+ try:
+ _ = dynsys.DynamicSystem(
+ time_interval_durations=time_step_durations,
+ A=B,
+ B=A,
+ C=C,
+ D=D,
+ integrate_outputs=integrate_y)
+ except ValueError:
+ number_errors += 1
+
+ assert number_errors == 1
+
+ # test C
+
+ number_errors = 0
+
+ try:
+ _ = dynsys.DynamicSystem(
+ time_interval_durations=time_step_durations,
+ A=C,
+ B=B,
+ C=A,
+ D=D,
+ integrate_outputs=integrate_y)
+ except ValueError:
+ number_errors += 1
+
+ assert number_errors == 1
+
+ # test D
+
+ number_errors = 0
+
+ try:
+ _ = dynsys.DynamicSystem(
+ time_interval_durations=time_step_durations,
+ A=D,
+ B=B,
+ C=C,
+ D=A,
+ integrate_outputs=integrate_y)
+ except ValueError:
+ number_errors += 1
+
+ assert number_errors == 1
+
+ # multiple matrices and time intervals, but more matrices than time steps
+
+ # *************************************************************************
+ # *************************************************************************
+
+ def test_unrecognised_matrix_formats(self):
+
+ integrate_outputs = False
+
+ # single matrix: matrix as lists of lists
+
+ time_step_durations = 8
+
+ A = 3
+
+ B = np.array([[3]])
+
+ C = np.array([[2]])
+
+ D = np.array([[8]])
+
+ # test A
+
+ number_errors = 0
+
+ try:
+ _ = dynsys.DynamicSystem(
+ time_interval_durations=time_step_durations,
+ A=A,
+ B=B,
+ C=C,
+ D=D,
+ integrate_outputs=integrate_outputs)
+ except TypeError:
+ number_errors += 1
+
+ assert number_errors == 1
+
+ # test B
+
+ number_errors = 0
+
+ try:
+ _ = dynsys.DynamicSystem(
+ time_interval_durations=time_step_durations,
+ A=B,
+ B=A,
+ C=C,
+ D=D,
+ integrate_outputs=integrate_outputs)
+ except TypeError:
+ number_errors += 1
+
+ assert number_errors == 1
+
+ # test C
+
+ number_errors = 0
+
+ try:
+ _ = dynsys.DynamicSystem(
+ time_interval_durations=time_step_durations,
+ A=C,
+ B=B,
+ C=A,
+ D=D,
+ integrate_outputs=integrate_outputs)
+ except TypeError:
+ number_errors += 1
+
+ assert number_errors == 1
+
+ # test D
+
+ number_errors = 0
+
+ try:
+ _ = dynsys.DynamicSystem(
+ time_interval_durations=time_step_durations,
+ A=D,
+ B=B,
+ C=C,
+ D=A,
+ integrate_outputs=integrate_outputs)
+ except TypeError:
+ number_errors += 1
+
+ assert number_errors == 1
+
+ # *********************************************************************
+
+ # multiple matrices: matrix as lists of lists
+
+ time_step_durations = [1, 1]
+
+ A = [[[3]], [[3]]]
+
+ B = [np.array([[3]]),np.array([[4]])]
+
+ C = [np.array([[2]]),np.array([[2]])]
+
+ D = [np.array([[8]]),np.array([[7]])]
+
+ # test A
+
+ number_errors = 0
+
+ try:
+ _ = dynsys.DynamicSystem(
+ time_interval_durations=time_step_durations,
+ A=A,
+ B=B,
+ C=C,
+ D=D,
+ integrate_outputs=integrate_outputs)
+ except TypeError:
+ number_errors += 1
+
+ assert number_errors == 1
+
+ # test B
+
+ number_errors = 0
+
+ try:
+ _ = dynsys.DynamicSystem(
+ time_interval_durations=time_step_durations,
+ A=B,
+ B=A,
+ C=C,
+ D=D,
+ integrate_outputs=integrate_outputs)
+ except TypeError:
+ number_errors += 1
+
+ assert number_errors == 1
+
+ # test C
+
+ number_errors = 0
+
+ try:
+ _ = dynsys.DynamicSystem(
+ time_interval_durations=time_step_durations,
+ A=C,
+ B=B,
+ C=A,
+ D=D,
+ integrate_outputs=integrate_outputs)
+ except TypeError:
+ number_errors += 1
+
+ assert number_errors == 1
+
+ # test D
+
+ number_errors = 0
+
+ try:
+ _ = dynsys.DynamicSystem(
+ time_interval_durations=time_step_durations,
+ A=D,
+ B=B,
+ C=C,
+ D=A,
+ integrate_outputs=integrate_outputs)
+ except TypeError:
+ number_errors += 1
+
+ assert number_errors == 1
+
+ # *********************************************************************
+ # *********************************************************************
+
+ # *************************************************************************
+ # *************************************************************************
+
+ def test_varying_matrix_sizes_no_integration(self):
+
+ integrate_outputs = False
+
+ time_step_durations = [1, 1]
+
+ A1 = np.array([[1]])
+ A2 = np.array([[1,3],[4,7]])
+
+ A = [A1, A2]
+ B = [np.array([[3]]),np.array([[4]])]
+ C = [np.array([[2]]),np.array([[2]])]
+ D = [np.array([[8]]),np.array([[7]])]
+
+ # test A
+ number_errors = 0
+ try:
+ _ = dynsys.DynamicSystem(
+ time_interval_durations=time_step_durations,
+ A=A,
+ B=B,
+ C=C,
+ D=D,
+ integrate_outputs=integrate_outputs)
+ except ValueError:
+ number_errors += 1
+ assert number_errors == 1
+
+ # test B
+ number_errors = 0
+ try:
+ _ = dynsys.DynamicSystem(
+ time_interval_durations=time_step_durations,
+ A=B,
+ B=A,
+ C=C,
+ D=D,
+ integrate_outputs=integrate_outputs)
+ except ValueError:
+ number_errors += 1
+ assert number_errors == 1
+
+ # test C
+ number_errors = 0
+ try:
+ _ = dynsys.DynamicSystem(
+ time_interval_durations=time_step_durations,
+ A=C,
+ B=B,
+ C=A,
+ D=D,
+ integrate_outputs=integrate_outputs)
+ except ValueError:
+ number_errors += 1
+ assert number_errors == 1
+
+ # test D
+ number_errors = 0
+ try:
+ _ = dynsys.DynamicSystem(
+ time_interval_durations=time_step_durations,
+ A=D,
+ B=B,
+ C=C,
+ D=A,
+ integrate_outputs=integrate_outputs)
+ except ValueError:
+ number_errors += 1
+ assert number_errors == 1
+
+ # *************************************************************************
+ # *************************************************************************
+
+ def test_varying_matrix_sizes_integration(self):
+
+ integrate_outputs = True
+
+ time_step_durations = [1, 1]
+
+ A1 = np.array([[1]])
+ A2 = np.array([[1,3],[4,7]])
+
+ A = [A1, A2]
+ B = [np.array([[3]]),np.array([[4]])]
+ C = [np.array([[2]]),np.array([[2]])]
+ D = [np.array([[8]]),np.array([[7]])]
+
+ # test A
+
+ number_errors = 0
+ try:
+ _ = dynsys.DynamicSystem(
+ time_interval_durations=time_step_durations,
+ A=A,
+ B=B,
+ C=C,
+ D=D,
+ integrate_outputs=integrate_outputs)
+ except ValueError:
+ number_errors += 1
+ assert number_errors == 1
+
+ # test B
+
+ number_errors = 0
+ try:
+ _ = dynsys.DynamicSystem(
+ time_interval_durations=time_step_durations,
+ A=B,
+ B=A,
+ C=C,
+ D=D,
+ integrate_outputs=integrate_outputs)
+ except ValueError:
+ number_errors += 1
+ assert number_errors == 1
+
+ # test C
+ number_errors = 0
+ try:
+ _ = dynsys.DynamicSystem(
+ time_interval_durations=time_step_durations,
+ A=C,
+ B=B,
+ C=A,
+ D=D,
+ integrate_outputs=integrate_outputs)
+ except ValueError:
+ number_errors += 1
+ assert number_errors == 1
+
+ # test D
+
+ number_errors = 0
+ try:
+ _ = dynsys.DynamicSystem(
+ time_interval_durations=time_step_durations,
+ A=D,
+ B=B,
+ C=C,
+ D=A,
+ integrate_outputs=integrate_outputs)
+ except ValueError:
+ number_errors += 1
+ assert number_errors == 1
+
+ # *************************************************************************
+ # *************************************************************************
+
+ def test_dynsys_stateless_multiout(self):
+
+ integrate_outputs = False
+ number_outputs = 2
+ method_dynsys_stateless_multiout(integrate_outputs, number_outputs)
+
+ integrate_outputs = False
+ number_outputs = 1
+ method_dynsys_stateless_multiout(integrate_outputs, number_outputs)
+
+ integrate_outputs = True
+ number_outputs = 2
+ method_dynsys_stateless_multiout(integrate_outputs, number_outputs)
+
+ integrate_outputs = True
+ number_outputs = 1
+ method_dynsys_stateless_multiout(integrate_outputs, number_outputs)
+
+ # integrate_outputs = False
+ # number_outputs = 0
+ # method_dynsys_stateless_multiout(integrate_outputs, number_outputs)
+
+ # integrate_outputs = True
+ # number_outputs = 0
+ # method_dynsys_stateless_multiout(integrate_outputs, number_outputs)
+
+ # *************************************************************************
+ # *************************************************************************
+
+ def test_dynsys_multiode_multiout(self):
+
+ integrate_outputs = False
+ number_outputs = 0
+ method_dynsys_multiode_multiout(integrate_outputs, number_outputs)
+
+ integrate_outputs = True
+ number_outputs = 0
+ method_dynsys_multiode_multiout(integrate_outputs, number_outputs)
+
+ integrate_outputs = False
+ number_outputs = 1
+ method_dynsys_multiode_multiout(integrate_outputs, number_outputs)
+
+ integrate_outputs = True
+ number_outputs = 1
+ method_dynsys_multiode_multiout(integrate_outputs, number_outputs)
+
+ integrate_outputs = False
+ number_outputs = 2
+ method_dynsys_multiode_multiout(integrate_outputs, number_outputs)
+
+ integrate_outputs = True
+ number_outputs = 2
+ method_dynsys_multiode_multiout(integrate_outputs, number_outputs)
+
+ # *************************************************************************
+ # *************************************************************************
+
+ def test_dynsys_singleode_multiout(self):
+
+ integrate_outputs = True
+ number_outputs = 2
+ method_dynsys_singleode_multiout(integrate_outputs, number_outputs)
+
+ integrate_outputs = True
+ number_outputs = 1
+ method_dynsys_singleode_multiout(integrate_outputs, number_outputs)
+
+ integrate_outputs = True
+ number_outputs = 0
+ method_dynsys_singleode_multiout(integrate_outputs, number_outputs)
+
+ integrate_outputs = False
+ number_outputs = 2
+ method_dynsys_singleode_multiout(integrate_outputs, number_outputs)
+
+ integrate_outputs = False
+ number_outputs = 1
+ method_dynsys_singleode_multiout(integrate_outputs, number_outputs)
+
+ integrate_outputs = False
+ number_outputs = 0
+ method_dynsys_singleode_multiout(integrate_outputs, number_outputs)
+
+ # *************************************************************************
+ # *************************************************************************
+
+def method_dynsys_stateless_multiout(integrate_outputs, number_outputs):
+
+ # *************************************************************************
+ # *************************************************************************
+
+ # time steps
+ list_regular_time_steps = [1, 1, 1, 1]
+ list_irregular_time_steps = [1, 1.5, 0.5, 1]
+ assert len(list_regular_time_steps) == len(list_irregular_time_steps)
+
+ # define the inputs
+ list_inputs = [-5, 50, 0.1, 1]
+
+ # *************************************************************************
+ # *************************************************************************
+
+ # generate time varying problem
+
+ list_A_matrices = []
+ list_B_matrices = []
+ list_C_matrices = []
+ list_D_matrices = []
+ for dt in list_irregular_time_steps:
+ # data
+ Aw, min_rel_heat = get_stateless_model_data(
+ relative_amplitude_variation=0.1)
+ # matrices
+ (A_matrix,
+ B_matrix,
+ C_matrix,
+ D_matrix) = stateless_model(Aw, min_rel_heat)
+ if number_outputs == 0:
+ D_matrix = None
+ elif number_outputs != None:
+ D_matrix = D_matrix[0:number_outputs,:]
+
+ list_A_matrices.append(A_matrix)
+ list_B_matrices.append(B_matrix)
+ list_C_matrices.append(C_matrix)
+ list_D_matrices.append(D_matrix)
+
+ if number_outputs == 0:
+ list_D_matrices = None
+
+ # *************************************************************************
+ # *************************************************************************
+
+ # generate time invariant problem
+
+ # data
+
+ Aw, min_rel_heat = get_stateless_model_data()
+
+ # matrices
+
+ (A_matrix,
+ B_matrix,
+ C_matrix,
+ D_matrix) = stateless_model(Aw, min_rel_heat)
+
+ if number_outputs == 0:
+ D_matrix = None
+ elif number_outputs != None:
+ D_matrix = D_matrix[0:number_outputs,:]
+
+ # *************************************************************************
+ # *************************************************************************
+
+ # time invariant model, regular time steps
+ x, y = example_dynsys_time_invariant(
+ list_regular_time_steps,
+ list_inputs,
+ integrate_outputs,
+ A_matrix,
+ B_matrix,
+ C_matrix,
+ D_matrix,
+ x0=None)
+
+ # TODO: check results
+
+ # *************************************************************************
+ # *************************************************************************
+
+ # time invariant model, irregular time steps
+ x, y = example_dynsys_time_invariant(
+ list_irregular_time_steps,
+ list_inputs,
+ integrate_outputs,
+ A_matrix,
+ B_matrix,
+ C_matrix,
+ D_matrix,
+ x0=None)
+
+ # TODO: check results
+
+# *****************************************************************************
+# *****************************************************************************
+
+def method_stateless_system_object(integrate_outputs: bool):
+
+ # regular time steps
+ time_step_durations = [1, 1, 1, 1]
+
+ # define a sequence of time steps
+ t = np.array([sum(time_step_durations[0:i])
+ for i in range(len(time_step_durations)+1)])
+
+ # define the inputs
+ list_inputs = [-5, 50, 0.1, 1]
+ U = np.array([[input_i for _ in range(len(time_step_durations))]
+ for input_i in list_inputs])
+
+ # get the model
+ Aw, min_rel_heat = get_stateless_model_data()
+ _, _, _, D = stateless_model(Aw, min_rel_heat)
+ ds = dynsys.StatelessSystem(
+ time_interval_durations=time_step_durations,
+ D=D,
+ integrate_outputs=integrate_outputs)
+
+ # define the initial conditions
+ X, Y = ds.simulate(U)
+ assert X == None
+ assert isinstance(Y, np.ndarray)
+ if integrate_outputs:
+ assert Y.shape == (2,len(t)-1) # two outputs
+ else:
+ assert Y.shape == (2,len(t)) # two outputs
+
+ # *********************************************************************
+
+ # regular time steps (as an int)
+ time_step_durations = [1]
+
+ # define a sequence of time steps
+ t = np.array([sum(time_step_durations[0:i])
+ for i in range(len(time_step_durations)+1)])
+
+ # define the inputs
+ list_inputs = [-5, 50, 0.1, 1]
+ U = np.array([[input_i for _ in range(len(time_step_durations))]
+ for input_i in list_inputs])
+
+ # get the model
+ Aw, min_rel_heat = get_stateless_model_data()
+ _, _, _, D = stateless_model(Aw, min_rel_heat)
+ ds = dynsys.StatelessSystem(
+ time_interval_durations=time_step_durations[0],
+ D=D,
+ integrate_outputs=integrate_outputs)
+
+ # define the initial conditions
+ X, Y = ds.simulate(U)
+ # X, Y = ds.simulate(U[:,:-1]) # to avoid having the second
+
+ assert X == None
+ assert isinstance(Y, np.ndarray)
+ if integrate_outputs:
+ assert Y.shape == (2,len(t)-1) # two outputs
+ else:
+ assert Y.shape == (2,len(t)) # two outputs
+
+ # *********************************************************************
+
+ # irregular time steps
+ time_step_durations = [1, 1.5, 0.5, 1]
+ # define a sequence of time steps
+ t = np.array([sum(time_step_durations[0:i])
+ for i in range(len(time_step_durations)+1)])
+ # define the inputs
+ list_inputs = [-5, 50, 0.1, 1]
+ U = np.array([[input_i for _ in range(len(time_step_durations))]
+ for input_i in list_inputs])
+ # get the model
+ Aw, min_rel_heat = get_stateless_model_data()
+ _, _, _, D = stateless_model(Aw, min_rel_heat)
+ ds = dynsys.StatelessSystem(
+ time_interval_durations=time_step_durations,
+ D=D,
+ integrate_outputs=integrate_outputs)
+ # define the initial conditions
+ X, Y = ds.simulate(U)
+ assert X == None
+ assert isinstance(Y, np.ndarray)
+ if integrate_outputs:
+ assert Y.shape == (2,len(t)-1) # two outputs
+ else:
+ assert Y.shape == (2,len(t)) # two outputs
+
+# *****************************************************************************
+# *****************************************************************************
+
+def method_dynsys_multiode_multiout(
+ integrate_outputs: bool,
+ number_outputs: int
+ ):
+
+ # *************************************************************************
+
+ # time steps
+ list_regular_time_steps = [1, 1, 1, 1]
+ list_irregular_time_steps = [1, 1.5, 0.5, 1]
+ assert len(list_regular_time_steps) == len(list_irregular_time_steps)
+
+ # define the inputs
+ list_inputs = [-5, 50, 0.1, 1]
+
+ # *************************************************************************
+
+ # generate time varying problem
+ list_A_matrices = []
+ list_B_matrices = []
+ list_C_matrices = []
+ list_D_matrices = []
+
+ for dt in list_irregular_time_steps:
+
+ # data
+ Ci, Ch, Ria, Rih, Aw, min_rel_heat, x0 = get_multi_ode_model_data(
+ relative_amplitude_variation=0.1)
+
+ # matrices
+ (A_matrix,
+ B_matrix,
+ C_matrix,
+ D_matrix) = two_node_model(Ci, Ch, Ria, Rih, Aw, min_rel_heat)
+
+ if number_outputs == 0:
+ C_matrix = None
+ D_matrix = None
+
+ elif number_outputs != None:
+ C_matrix = C_matrix[0:number_outputs,:]
+ D_matrix = D_matrix[0:number_outputs,:]
+
+ list_A_matrices.append(A_matrix)
+ list_B_matrices.append(B_matrix)
+ list_C_matrices.append(C_matrix)
+ list_D_matrices.append(D_matrix)
+
+ if number_outputs == 0:
+ list_C_matrices = None
+ list_D_matrices = None
+
+ # *************************************************************************
+
+ # generate time invariant problem
+
+ # data
+ Ci, Ch, Ria, Rih, Aw, min_rel_heat, x0 = get_multi_ode_model_data()
+
+ # matrices
+ (A_matrix,
+ B_matrix,
+ C_matrix,
+ D_matrix) = two_node_model(Ci, Ch, Ria, Rih, Aw, min_rel_heat)
+
+ if number_outputs == 0:
+ C_matrix = None
+ D_matrix = None
+ elif number_outputs != None:
+ C_matrix = C_matrix[0:number_outputs,:]
+ D_matrix = D_matrix[0:number_outputs,:]
+
+ # *************************************************************************
+
+ # time invariant model, regular time steps
+
+ x_scipy, y_scipy = example_scipy_time_invariant_regular_steps(
+ list_regular_time_steps,
+ list_inputs,
+ integrate_outputs,
+ A_matrix,
+ B_matrix,
+ C_matrix,
+ D_matrix,
+ x0)
+
+ x, y = example_dynsys_time_invariant(
+ list_regular_time_steps,
+ list_inputs,
+ integrate_outputs,
+ A_matrix,
+ B_matrix,
+ C_matrix,
+ D_matrix,
+ x0)
+
+ for (x_scipy_i, x_i) in zip(x_scipy, x):
+ np.testing.assert_allclose(x_i, x_scipy_i)
+
+ if number_outputs != 0:
+ for (y_scipy_i, y_i) in zip(y_scipy, y):
+ np.testing.assert_allclose(y_i, y_scipy_i)
+
+ # *************************************************************************
+
+ # time invariant model, irregular time steps
+
+ x_scipy, y_scipy = example_scipy_time_invariant_irregular_steps(
+ list_irregular_time_steps,
+ list_inputs,
+ integrate_outputs,
+ A_matrix,
+ B_matrix,
+ C_matrix,
+ D_matrix,
+ x0)
+
+ x, y = example_dynsys_time_invariant(
+ list_irregular_time_steps,
+ list_inputs,
+ integrate_outputs,
+ A_matrix,
+ B_matrix,
+ C_matrix,
+ D_matrix,
+ x0)
+
+ for (x_scipy_i, x_i) in zip(x_scipy, x):
+ np.testing.assert_allclose(x_i, x_scipy_i)
+
+ if number_outputs != 0:
+ for (y_scipy_i, y_i) in zip(y_scipy, y):
+ np.testing.assert_allclose(y_i, y_scipy_i)
+
+ # *************************************************************************
+
+ # time-varying model, regular time steps
+
+ x_scipy, y_scipy = example_scipy_time_varying(
+ list_regular_time_steps,
+ list_inputs,
+ integrate_outputs,
+ list_A_matrices,
+ list_B_matrices,
+ list_C_matrices,
+ list_D_matrices,
+ x0)
+
+ x, y = example_dynsys_time_varying(
+ list_regular_time_steps,
+ list_inputs,
+ integrate_outputs,
+ list_A_matrices,
+ list_B_matrices,
+ list_C_matrices,
+ list_D_matrices,
+ x0)
+
+ for (x_scipy_i, x_i) in zip(x_scipy, x):
+ np.testing.assert_allclose(x_i, x_scipy_i)
+
+ if number_outputs != 0:
+ for (y_scipy_i, y_i) in zip(y_scipy, y):
+ np.testing.assert_allclose(y_i, y_scipy_i)
+
+ # *************************************************************************
+
+ # time-varying model, irregular time steps
+
+ x_scipy, y_scipy = example_scipy_time_varying(
+ list_irregular_time_steps,
+ list_inputs,
+ integrate_outputs,
+ list_A_matrices,
+ list_B_matrices,
+ list_C_matrices,
+ list_D_matrices,
+ x0)
+
+ x, y = example_dynsys_time_varying(
+ list_irregular_time_steps,
+ list_inputs,
+ integrate_outputs,
+ list_A_matrices,
+ list_B_matrices,
+ list_C_matrices,
+ list_D_matrices,
+ x0)
+
+ for (x_scipy_i, x_i) in zip(x_scipy, x):
+ np.testing.assert_allclose(x_i, x_scipy_i)
+
+ if number_outputs != 0:
+ for (y_scipy_i, y_i) in zip(y_scipy, y):
+ np.testing.assert_allclose(y_i, y_scipy_i)
+
+# *****************************************************************************
+# *****************************************************************************
+
+# test single ODE, multi-output dynamic systems while integrating outputs
+
+# examples_dynsys_singleode_multiout(True, print_outputs, seed_number)
+
+# # test single ODE, multi-output dynamic systems without integrating outputs
+
+# examples_dynsys_singleode_multiout(False, print_outputs, seed_number)
+
+# # test single-ODE, single-output dynamic systems while integrating outputs
+
+# examples_dynsys_singleode_multiout(True, print_outputs, seed_number, 1)
+
+# # test single-ODE, single-output dynamic systems without integrating outputs
+
+# examples_dynsys_singleode_multiout(False, print_outputs, seed_number, 1)
+
+# # *************************************************************************
+
+# # outputless
+
+# # test single-ODE, outputless dynamic systems while integrating outputs
+
+# examples_dynsys_singleode_multiout(True, print_outputs, seed_number, 0)
+
+# # test single-ODE, outputless dynamic systems without integrating outputs
+
+# examples_dynsys_singleode_multiout(False, print_outputs, seed_number, 0)
+
+def method_dynsys_singleode_multiout(integrate_outputs: bool,
+ number_outputs: int = 2):
+
+ # *************************************************************************
+
+ # time steps
+ list_regular_time_steps = [1, 1, 1, 1]
+ list_irregular_time_steps = [1, 1.5, 0.5, 1]
+ assert len(list_regular_time_steps) == len(list_irregular_time_steps)
+
+ # define the inputs
+ list_inputs = [-5, 50, 0.1, 1]
+
+ # *************************************************************************
+
+ # generate time varying problem
+ list_A_matrices = []
+ list_B_matrices = []
+ list_C_matrices = []
+ list_D_matrices = []
+
+ for dt in list_irregular_time_steps:
+
+ # data
+
+ Ci, Ria, Aw, min_rel_heat, x0 = get_single_ode_model_data(
+ relative_amplitude_variation=0.1)
+
+ # matrices
+
+ (A_matrix,
+ B_matrix,
+ C_matrix,
+ D_matrix) = single_node_model(Ci, Ria, Aw, min_rel_heat)
+
+ if number_outputs == 0:
+ C_matrix = None
+ D_matrix = None
+ elif number_outputs != None:
+ C_matrix = C_matrix[0:number_outputs,:]
+ D_matrix = D_matrix[0:number_outputs,:]
+
+ list_A_matrices.append(A_matrix)
+ list_B_matrices.append(B_matrix)
+ list_C_matrices.append(C_matrix)
+ list_D_matrices.append(D_matrix)
+
+ if number_outputs == 0:
+ list_C_matrices = None
+ list_D_matrices = None
+
+ # *************************************************************************
+
+ # generate time invariant problem
+
+ # data
+
+ Ci, Ria, Aw, min_rel_heat, x0 = get_single_ode_model_data()
+
+ # matrices
+
+ (A_matrix,
+ B_matrix,
+ C_matrix,
+ D_matrix) = single_node_model(Ci, Ria, Aw, min_rel_heat)
+
+ if number_outputs == 0:
+ C_matrix = None
+ D_matrix = None
+ elif number_outputs != None:
+ C_matrix = C_matrix[0:number_outputs,:]
+ D_matrix = D_matrix[0:number_outputs,:]
+
+ # *************************************************************************
+
+ # time invariant model, regular time steps
+
+ x_scipy, y_scipy = example_scipy_time_invariant_regular_steps(
+ list_regular_time_steps,
+ list_inputs,
+ integrate_outputs,
+ A_matrix,
+ B_matrix,
+ C_matrix,
+ D_matrix,
+ x0)
+
+ x, y = example_dynsys_time_invariant(
+ list_regular_time_steps,
+ list_inputs,
+ integrate_outputs,
+ A_matrix,
+ B_matrix,
+ C_matrix,
+ D_matrix,
+ x0)
+
+ for (x_scipy_i, x_i) in zip(x_scipy, x):
+ np.testing.assert_allclose(x_i, x_scipy_i)
+
+ if number_outputs != 0:
+ for (y_scipy_i, y_i) in zip(y_scipy, y):
+ np.testing.assert_allclose(y_i, y_scipy_i)
+
+ # *************************************************************************
+
+ # time invariant model, irregular time steps
+
+ x_scipy, y_scipy = example_scipy_time_invariant_irregular_steps(
+ list_irregular_time_steps,
+ list_inputs,
+ integrate_outputs,
+ A_matrix,
+ B_matrix,
+ C_matrix,
+ D_matrix,
+ x0)
+
+ x, y = example_dynsys_time_invariant(
+ list_irregular_time_steps,
+ list_inputs,
+ integrate_outputs,
+ A_matrix,
+ B_matrix,
+ C_matrix,
+ D_matrix,
+ x0)
+
+ for (x_scipy_i, x_i) in zip(x_scipy, x):
+ np.testing.assert_allclose(x_i, x_scipy_i)
+
+ if number_outputs != 0:
+ for (y_scipy_i, y_i) in zip(y_scipy, y):
+ np.testing.assert_allclose(y_i, y_scipy_i)
+
+ # *************************************************************************
+
+ # time-varying model, regular time steps
+
+ x_scipy, y_scipy = example_scipy_time_varying(
+ list_regular_time_steps,
+ list_inputs,
+ integrate_outputs,
+ list_A_matrices,
+ list_B_matrices,
+ list_C_matrices,
+ list_D_matrices,
+ x0)
+
+ x, y = example_dynsys_time_varying(
+ list_regular_time_steps,
+ list_inputs,
+ integrate_outputs,
+ list_A_matrices,
+ list_B_matrices,
+ list_C_matrices,
+ list_D_matrices,
+ x0)
+
+ for (x_scipy_i, x_i) in zip(x_scipy, x):
+ np.testing.assert_allclose(x_i, x_scipy_i)
+
+ if number_outputs != 0:
+ for (y_scipy_i, y_i) in zip(y_scipy, y):
+ np.testing.assert_allclose(y_i, y_scipy_i)
+
+ # *************************************************************************
+
+ # time-varying model, irregular time steps
+
+ x_scipy, y_scipy = example_scipy_time_varying(
+ list_irregular_time_steps,
+ list_inputs,
+ integrate_outputs,
+ list_A_matrices,
+ list_B_matrices,
+ list_C_matrices,
+ list_D_matrices,
+ x0)
+
+ x, y = example_dynsys_time_varying(
+ list_irregular_time_steps,
+ list_inputs,
+ integrate_outputs,
+ list_A_matrices,
+ list_B_matrices,
+ list_C_matrices,
+ list_D_matrices,
+ x0)
+
+ for (x_scipy_i, x_i) in zip(x_scipy, x):
+ np.testing.assert_allclose(x_i, x_scipy_i)
+
+ if number_outputs != 0:
+ for (y_scipy_i, y_i) in zip(y_scipy, y):
+ np.testing.assert_allclose(y_i, y_scipy_i)
+
+# *****************************************************************************
+# *****************************************************************************
+
+def example_scipy_time_invariant_regular_steps(
+ time_step_durations: list,
+ inputs_list: list,
+ integrate_outputs: bool,
+ A_matrix,
+ B_matrix,
+ C_matrix,
+ D_matrix,
+ x0):
+
+ # define a sequence of time steps
+
+ t = np.array([sum(time_step_durations[0:i])
+ for i in range(len(time_step_durations)+1)])
+
+ # define the inputs
+
+ u = np.array([inputs_list for dt in t])
+
+ if type(C_matrix) == type(None) and type(D_matrix) == type(None):
+
+ ss = StateSpace(A_matrix,
+ B_matrix,
+ np.zeros((1,A_matrix.shape[1])),
+ np.zeros((1,B_matrix.shape[1])))
+
+ else:
+
+ ss = StateSpace(A_matrix, B_matrix, C_matrix, D_matrix)
+
+ # simulate the system response
+
+ tout, yout, xout = lsim(ss, U=u, T=t, X0=x0)
+
+ # convert to matrix format if there is only one dimension
+
+ if len(xout.shape) == 1:
+
+ xout = np.array([xout]).T
+
+ if len(yout.shape) == 1:
+
+ yout = np.array([yout]).T
+
+ if integrate_outputs:
+
+ # yout's shape should be: (number of time steps, 2)
+
+ # ignore the first instant
+
+ yout = yout[1:,:]
+
+ yout_m, yout_n = yout.shape
+
+ # multiply each output by the respective time step
+
+ yout = np.array(
+ [[yout[m,n]*time_step_durations[n] for n in range(yout_n)]
+ for m in range(yout_m)]
+ )
+
+ # note: the code above should not work when C != 0
+
+ return xout.T, yout.T
+
+# *****************************************************************************
+# *****************************************************************************
+
+def example_scipy_time_invariant_irregular_steps(
+ time_step_durations: list,
+ inputs_list: list,
+ integrate_outputs: bool,
+ A_matrix,
+ B_matrix,
+ C_matrix,
+ D_matrix,
+ x0):
+
+ # define a sequence of time steps
+
+ t = np.array([sum(time_step_durations[0:i])
+ for i in range(len(time_step_durations)+1)])
+
+ # define the inputs
+
+ u = np.array([inputs_list for _ in t])
+
+ # get the model
+
+ if type(C_matrix) == type(None) and type(D_matrix) == type(None):
+
+ ss = StateSpace(A_matrix,
+ B_matrix,
+ np.zeros((1,A_matrix.shape[1])),
+ np.zeros((1,B_matrix.shape[1])))
+
+ yout = np.zeros((1, len(time_step_durations)+1))
+
+ else:
+
+ ss = StateSpace(A_matrix, B_matrix, C_matrix, D_matrix)
+
+ yout = np.zeros((C_matrix.shape[0], len(time_step_durations)+1))
+
+ # declare output variables
+
+ xout = np.zeros((A_matrix.shape[0], len(time_step_durations)+1))
+
+ #yout = np.zeros((C_matrix.shape[0], len(time_step_durations)+1))
+
+ xout[:,0] = x0
+
+ # initial output
+
+ if not integrate_outputs:
+
+ # do not ignore the first instant
+
+ yout[:,0] = np.dot(ss.D, u[0,:])
+
+ # else:
+
+ # yout[:,0] = np.array([0, 0])
+
+ # for each time step
+
+ for i in range(len(time_step_durations)):
+
+ # time step duration
+
+ dt = time_step_durations[i]
+
+ ssd = ss.to_discrete(dt=dt)
+
+ # simulate the system response
+
+ _, ytemp, xtemp = dlsim(ssd,
+ u=np.array([u[i,:],u[i+1,:]]),
+ t=np.array([0,dt]),
+ x0=xout[:,i])
+
+ xout[:,i+1] = xtemp[1,:] #np.array([0,0])
+
+ yout[:,i+1] = ytemp[1,:] #np.array([0,0])
+
+ if integrate_outputs:
+
+ # yout's shape should be: (2,number of time steps)
+
+ # ignore the first instant
+
+ yout = yout[:,1:]
+
+ yout_m, yout_n = yout.shape
+
+ # multiply each output by the respective time step
+
+ for m in range(yout_m):
+
+ for n in range(yout_n):
+
+ if n == 0:
+
+ continue # ignore the first time interval
+
+ yout[m,n] = yout[m,n]*time_step_durations[n]
+
+ return xout, yout
+
+# *****************************************************************************
+# *****************************************************************************
+
+def example_scipy_time_varying(time_step_durations: list,
+ inputs_list: list,
+ integrate_outputs: bool,
+ A_matrix,
+ B_matrix,
+ C_matrix,
+ D_matrix,
+ x0):
+
+ # define a sequence of time steps
+
+ t = np.array([sum(time_step_durations[0:i])
+ for i in range(len(time_step_durations)+1)])
+
+ # define the inputs
+
+ u = np.array([inputs_list for dt in t])
+
+ # define the initial conditions
+
+ xout = np.zeros((A_matrix[0].shape[0], len(time_step_durations)+1))
+
+ xout[:,0] = x0
+
+ if type(C_matrix) == type(None) and type(D_matrix) == type(None):
+
+ yout = np.zeros((1, len(time_step_durations)+1))
+
+ else:
+
+ yout = np.zeros((C_matrix[0].shape[0], len(time_step_durations)+1))
+
+
+ for i, dt in enumerate(time_step_durations):
+
+ # get the model
+
+ if type(C_matrix) == type(None) and type(D_matrix) == type(None):
+
+ ss = StateSpace(A_matrix[i],
+ B_matrix[i],
+ np.zeros((1,A_matrix[i].shape[1])),
+ np.zeros((1,B_matrix[i].shape[1])))
+
+ else:
+
+ ss = StateSpace(A_matrix[i], B_matrix[i], C_matrix[i], D_matrix[i])
+
+ if i == 0 and not integrate_outputs:
+
+ # compute the initial output
+
+ yout[:,0] = np.dot(ss.D, u[0,:])
+
+ ssd = ss.to_discrete(dt=dt)
+
+ # simulate the system response
+
+ if i == 0:
+
+ _, ytemp, xtemp = dlsim(ssd,
+ u=np.array([u[i],u[i]]),
+ t=np.array([0,dt]),
+ x0=x0)
+
+ else:
+
+ _, ytemp, xtemp = dlsim(ssd,
+ u=np.array([u[i],u[i]]),
+ t=np.array([0,dt]),
+ x0=xtemp[1,:])
+
+ # assert that the sizes match
+
+ xout[:,i+1] = xtemp[1,:]
+
+ yout[:,i+1] = ytemp[1,:]
+
+ if integrate_outputs:
+
+ # yout's shape should be: (2,number of time steps)
+
+ # ignore the first instant
+
+ yout = yout[:,1:]
+
+ yout_m, yout_n = yout.shape
+
+ # multiply each output by the respective time step
+
+ for m in range(yout_m):
+
+ for n in range(yout_n):
+
+ if n == 0:
+
+ continue # ignore the first time interval
+
+ yout[m,n] = yout[m,n]*time_step_durations[n]
+
+ return xout, yout
+
+# *****************************************************************************
+# *****************************************************************************
+
+def example_dynsys_time_invariant(time_step_durations: list,
+ inputs_list: list,
+ integrate_outputs: bool,
+ A_matrix,
+ B_matrix,
+ C_matrix,
+ D_matrix,
+ x0):
+
+ # define a sequence of time steps
+
+ t = np.array([sum(time_step_durations[0:i])
+ for i in range(len(time_step_durations)+1)])
+
+ # define the inputs
+
+ # U = np.array([[input_i for dt in t]
+ # for input_i in inputs_list])
+
+ U = np.array([[input_i for _ in range(len(time_step_durations))]
+ for input_i in inputs_list])
+
+ # get the model
+
+ ds = dynsys.DynamicSystem(
+ time_interval_durations=time_step_durations,
+ A=A_matrix,
+ B=B_matrix,
+ C=C_matrix,
+ D=D_matrix,
+ integrate_outputs=integrate_outputs)
+
+ # define the initial conditions
+
+ X, Y = ds.simulate(U, x0)
+
+ return X, Y
+
+# *****************************************************************************
+# *****************************************************************************
+
+def example_dynsys_time_varying(time_step_durations: list,
+ inputs_list: list,
+ integrate_outputs: bool,
+ A_matrix,
+ B_matrix,
+ C_matrix,
+ D_matrix,
+ x0):
+
+ # define a sequence of time steps
+
+ t = np.array([sum(time_step_durations[0:i])
+ for i in range(len(time_step_durations)+1)])
+
+ # define the inputs
+
+ # U = np.array([[input_i for dt in t]
+ # for input_i in inputs_list])
+
+ U = np.array([[input_i for _ in range(len(time_step_durations))]
+ for input_i in inputs_list])
+
+ ds = dynsys.DynamicSystem(
+ time_interval_durations=time_step_durations,
+ A=A_matrix,
+ B=B_matrix,
+ C=C_matrix,
+ D=D_matrix,
+ integrate_outputs=integrate_outputs)
+
+ # define the initial conditions
+
+ X, Y = ds.simulate(U, x0)
+
+ return X, Y
+
+# *****************************************************************************
+# *****************************************************************************
+
+def get_stateless_model_data(relative_amplitude_variation: float = 0.0):
+
+ mrh_deviation = random.random() - 0.5
+
+ Aw = 6.22 # original: 6.22 m2
+
+ min_rel_heat = 0.2*(1+relative_amplitude_variation*mrh_deviation)
+
+ return Aw, min_rel_heat
+
+# *****************************************************************************
+# *****************************************************************************
+
+def get_single_ode_model_data(relative_amplitude_variation: float = 0.0):
+
+ # define how the coefficients change
+
+ Ria_deviation = random.random() - 0.5
+
+ # define the (A, B, C and D) matrices
+ # A: n*n
+ # B: n*m
+ # C: r*n
+ # D: r*m
+
+ Ci = 1.360*3600000
+ Ria = (
+ (1+relative_amplitude_variation*Ria_deviation)*5.31/3600000
+ )
+ Aw = 6.22
+
+ min_rel_heat = 0.2
+
+ x0 = np.array([20])
+
+ return Ci, Ria, Aw, min_rel_heat, x0
+
+# *****************************************************************************
+# *****************************************************************************
+
+def get_multi_ode_model_data(relative_amplitude_variation: float = 0.0):
+
+ # define how the coefficients change
+
+ Rih_deviation = random.random() - 0.5
+
+ Ria_deviation = random.random() - 0.5
+
+ # define the (A, B, C and D) matrices
+ # A: n*n
+ # B: n*m
+ # C: r*n
+ # D: r*m
+
+ # from Bacher and Madsen (2011): model TiTh
+
+ Ci = 1.360*3600000 # original: 1.36 kWh/ºC
+ Ch = 0.309*3600000 # original: 0.309 kWh/ºC
+ Ria = (
+ (1+relative_amplitude_variation*Ria_deviation)*5.31/3600000
+ ) # original: 5.31 ºC/kWh
+ Rih = (
+ (1+relative_amplitude_variation*Rih_deviation)*0.639/3600000
+ ) # original: 0.639 ºC/kWh
+ Aw = 6.22 # original: 6.22 m2
+
+ min_rel_heat = 0.2
+
+ x0 = np.array([20, 20])
+
+ return Ci, Ch, Ria, Rih, Aw, min_rel_heat, x0
+
+# *****************************************************************************
+# *****************************************************************************
+
+def stateless_model(Aw, min_rel_heat):
+
+ # inputs: Ta, phi_s, phi_h above the minimum, phi_h status
+ # outputs: solar irradiance, heat
+
+ d = np.array([[0, Aw, 0, 0],
+ [0, 0, (1-min_rel_heat), min_rel_heat]])
+
+ return None, None, None, d
+
+# *****************************************************************************
+# *****************************************************************************
+
+def single_node_model(Ci, Ria, Aw, min_rel_heat):
+
+ # states: Ti and Th
+ # inputs: Ta, phi_s, phi_h above the minimum, phi_h status
+ # outputs: solar irradiance, heat
+
+ a = np.array([[-1/(Ria*Ci)]])
+ b = np.array([[1/(Ci*Ria), Aw/Ci, (1-min_rel_heat)/Ci, min_rel_heat/Ci]])
+ c = np.array([[0],[0]])
+ d = np.array([[0, Aw, 0, 0],
+ [0, 0, (1-min_rel_heat), min_rel_heat]])
+
+ return a, b, c, d
+
+# *****************************************************************************
+# *****************************************************************************
+
+def two_node_model(Ci, Ch, Ria, Rih, Aw, min_rel_heat):
+
+ # states: Ti and Th
+ # inputs: Ta, phi_s, phi_h above the minimum, phi_h status
+ # outputs: solar irradiance, heat
+
+ a = np.array([[-(1/Rih+1/Ria)/Ci, 1/(Ci*Rih)],
+ [1/(Ch*Rih), -1/(Ch*Rih)]])
+ b = np.array([[1/(Ci*Ria), Aw/Ci, 0, 0],
+ [0, 0, (1-min_rel_heat)/Ch, min_rel_heat/Ch]])
+ c = np.array([[0, 0],
+ [0, 0]])
+ d = np.array([[0, Aw, 0, 0],
+ [0, 0, (1-min_rel_heat), min_rel_heat]])
+
+ return a, b, c, d
+
+# *****************************************************************************
+# *****************************************************************************
\ No newline at end of file
diff --git a/tests/test_esipp_network.py b/tests/test_esipp_network.py
new file mode 100644
index 0000000..0959c95
--- /dev/null
+++ b/tests/test_esipp_network.py
@@ -0,0 +1,2061 @@
+# imports
+
+# standard
+
+import random
+
+from networkx import binomial_tree, MultiDiGraph
+
+# local
+
+from src.topupopt.problems.esipp.network import Arcs, Network
+
+from src.topupopt.problems.esipp.network import ArcsWithoutLosses
+
+from src.topupopt.problems.esipp.network import ArcsWithoutProportionalLosses
+
+from src.topupopt.problems.esipp.network import ArcsWithoutStaticLosses
+
+from src.topupopt.problems.esipp.resource import ResourcePrice
+
+# *****************************************************************************
+# *****************************************************************************
+
+# TODO: add test for directed arcs between import and export nodes with static losses
+
+# TODO: add test for undirected arcs involving import and export nodes
+
+class TestNetwork:
+
+ # *************************************************************************
+ # *************************************************************************
+
+ def test_tree_topology(self):
+
+ # create a network object with a tree topology
+
+ tree_network = binomial_tree(3, create_using=MultiDiGraph)
+
+ network = Network(tree_network)
+
+ for edge_key in network.edges(keys=True):
+
+ arc = ArcsWithoutLosses(
+ name=str(edge_key),
+ capacity=[5, 10],
+ minimum_cost=[3, 6],
+ specific_capacity_cost=0,
+ capacity_is_instantaneous=False
+ )
+
+ network.add_edge(
+ *edge_key,
+ **{Network.KEY_ARC_TECH: arc}
+ )
+
+ # assert that it does not have a tree topology
+
+ assert not network.has_tree_topology()
+
+ # select all the nodes
+
+ for edge_key in network.edges(keys=True):
+
+ network.edges[edge_key][
+ Network.KEY_ARC_TECH].options_selected[0] = True
+
+ # assert that it has a tree topology
+
+ assert network.has_tree_topology()
+
+ # *************************************************************************
+ # *************************************************************************
+
+ def test_arc_technologies_static_losses(self):
+
+ # *********************************************************************
+ # *********************************************************************
+
+ number_time_intervals = 3
+ number_scenarios = 2
+ number_options = 4
+
+ efficiency_dict = {
+ (q,k): 0.95
+ for q in range(number_scenarios)
+ for k in range(number_time_intervals)
+ }
+
+ static_loss_dict = {
+ (h,q,k): 1
+ for h in range(number_options)
+ for q in range(number_scenarios)
+ for k in range(number_time_intervals)
+ }
+
+ for capacity_is_instantaneous in (True, False):
+
+ arc_tech = Arcs(
+ name='any',
+ efficiency=efficiency_dict,
+ efficiency_reverse=None,
+ capacity=tuple(1+o for o in range(number_options)),
+ minimum_cost=tuple(1+o for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ static_loss=static_loss_dict,
+ validate=True
+ )
+
+ assert arc_tech.has_proportional_losses()
+
+ assert arc_tech.has_static_losses()
+
+ assert not arc_tech.is_infinite_capacity()
+
+ assert not arc_tech.has_been_selected()
+
+ assert arc_tech.is_isotropic(reverse_none_means_isotropic=True)
+
+ assert not arc_tech.is_isotropic(reverse_none_means_isotropic=False)
+
+ # isotropic
+
+ arc_tech = Arcs(
+ name='any',
+ efficiency=None,
+ efficiency_reverse=None,
+ capacity=tuple(1+o for o in range(number_options)),
+ minimum_cost=tuple(1+o for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ static_loss=static_loss_dict,
+ validate=True
+ )
+
+ assert not arc_tech.has_proportional_losses()
+
+ assert arc_tech.has_static_losses()
+
+ assert not arc_tech.is_infinite_capacity()
+
+ assert not arc_tech.has_been_selected()
+
+ assert arc_tech.is_isotropic(reverse_none_means_isotropic=True)
+
+ assert arc_tech.is_isotropic(reverse_none_means_isotropic=False)
+
+ # create arc technology with only one option
+
+ arc_tech = Arcs(
+ name='any',
+ efficiency=efficiency_dict,
+ efficiency_reverse=None,
+ capacity=(1,),
+ minimum_cost=(1,),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ static_loss={
+ (0,q,k): 1
+ #for h in range(number_options)
+ for q in range(number_scenarios)
+ for k in range(number_time_intervals)
+ },
+ validate=True
+ )
+
+ assert arc_tech.has_proportional_losses()
+
+ assert arc_tech.has_static_losses()
+
+ assert not arc_tech.is_infinite_capacity()
+
+ assert not arc_tech.has_been_selected()
+
+ assert arc_tech.is_isotropic(reverse_none_means_isotropic=True)
+
+ assert not arc_tech.is_isotropic(reverse_none_means_isotropic=False)
+
+ # create arc technology for one time interval
+
+ arc_tech = Arcs(
+ name='any',
+ efficiency={
+ (q, 0): 0.5
+ for q in range(number_scenarios)
+ #for k in range(number_time_intervals)
+ },
+ efficiency_reverse=None,
+ capacity=tuple(1+o for o in range(number_options)),
+ minimum_cost=tuple(1+o for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ static_loss={
+ (h,q,0): 1
+ for h in range(number_options)
+ for q in range(number_scenarios)
+ #for k in range(number_time_intervals)
+ },
+ validate=True
+ )
+
+ assert arc_tech.has_proportional_losses()
+
+ assert arc_tech.has_static_losses()
+
+ assert not arc_tech.is_infinite_capacity()
+
+ assert not arc_tech.has_been_selected()
+
+ assert arc_tech.is_isotropic(reverse_none_means_isotropic=True)
+
+ assert not arc_tech.is_isotropic(reverse_none_means_isotropic=False)
+
+ # *********************************************************************
+
+ # TypeError: The static losses should be given as a dict or None.
+
+ error_triggered = False
+ try:
+ _ = Arcs(
+ name='any',
+ efficiency=None,
+ efficiency_reverse=None,
+ capacity=tuple(1+o for o in range(number_options)),
+ minimum_cost=tuple(1+o for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ static_loss=tuple(
+ [k for k in range(number_time_intervals)]
+ for o in range(number_options)),
+ validate=True
+ )
+ except TypeError:
+ error_triggered = True
+ assert error_triggered
+
+ # ValueError('The static losses should be specified for each arc
+ # option.')
+
+ error_triggered = False
+ try:
+ _ = Arcs(
+ name='any',
+ efficiency=None,
+ efficiency_reverse=None,
+ capacity=tuple(1+o for o in range(number_options)),
+ minimum_cost=tuple(1+o for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ static_loss={
+ (h, q,): 1
+ for h in range(number_options)
+ for q in range(number_scenarios)
+ },
+ validate=True
+ )
+ except ValueError:
+ error_triggered = True
+ assert error_triggered
+
+ # TypeError('The static losses must be specified via a list of lists.')
+
+ error_triggered = False
+ try:
+ _ = Arcs(
+ name='any',
+ efficiency=None,
+ efficiency_reverse=None,
+ capacity=tuple(1+o for o in range(number_options)),
+ minimum_cost=tuple(1+o for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ static_loss=[
+ tuple(k for k in range(number_time_intervals))
+ for o in range(number_options)],
+ validate=True
+ )
+ except TypeError:
+ error_triggered = True
+ assert error_triggered
+
+ # ValueError('The static loss values are inconsistent with the number '
+ # 'of options, scenarios and intervals.')
+
+ error_triggered = False
+ try:
+ arc_tech = Arcs(
+ name='any',
+ efficiency=None,
+ efficiency_reverse=None,
+ capacity=tuple(1+o for o in range(number_options)),
+ minimum_cost=tuple(1+o for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ static_loss={
+ (h,q,k): 1
+ for h in range(number_options)
+ for q in range(number_scenarios)
+ for k in range(number_time_intervals-1)
+ },
+ validate=True
+ )
+
+ arc_tech.validate_sizes(number_options=number_options,
+ number_scenarios=number_scenarios,
+ number_intervals=[
+ number_time_intervals
+ for _ in range(number_scenarios)])
+ except ValueError:
+ error_triggered = True
+ assert error_triggered
+
+ # TypeError('The static losses were not provided as numbers.')
+
+ error_triggered = False
+ try:
+ _ = Arcs(
+ name='any',
+ efficiency=None,
+ efficiency_reverse=None,
+ capacity=tuple(1+o for o in range(number_options)),
+ minimum_cost=tuple(1+o for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ static_loss={
+ (h,q,k): str(3.54)
+ for h in range(number_options)
+ for q in range(number_scenarios)
+ for k in range(number_time_intervals)
+ },
+ validate=True
+ )
+ except TypeError:
+ error_triggered = True
+ assert error_triggered
+
+ # ValueError('The static losses must be positive or zero.')
+
+ error_triggered = False
+ try:
+ _ = Arcs(
+ name='any',
+ efficiency=None,
+ efficiency_reverse=None,
+ capacity=tuple(1+o for o in range(number_options)),
+ minimum_cost=tuple(1+o for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ static_loss={
+ (h,q,k): -random.randint(0, 1)*random.random()
+ for h in range(number_options)
+ for q in range(number_scenarios)
+ for k in range(number_time_intervals)
+ },
+ validate=True
+ )
+ except ValueError:
+ error_triggered = True
+ assert error_triggered
+
+ # TypeError: The static loss dict keys must be tuples
+
+ error_triggered = False
+ try:
+ _ = Arcs(
+ name='hey',
+ efficiency=None,
+ efficiency_reverse=None,
+ static_loss={k:1 for k in range(number_time_intervals)},
+ capacity=tuple(1+o for o in range(number_options)),
+ minimum_cost=tuple(1+o for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ validate=True)
+ except TypeError:
+ error_triggered = True
+ assert error_triggered
+
+ #ValueError( 'The static loss dict keys must be tuples of size 3.')
+
+ error_triggered = False
+ try:
+ _ = Arcs(
+ name='hey',
+ efficiency=None,
+ efficiency_reverse=None,
+ static_loss={(k,3): 1 for k in range(number_time_intervals)},
+ capacity=tuple(1+o for o in range(number_options)),
+ minimum_cost=tuple(1+o for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ validate=True)
+ except ValueError:
+ error_triggered = True
+ assert error_triggered
+
+ # TypeError(The staticl osses should be given as a dict or None.')
+
+ error_triggered = False
+ try:
+ _ = Arcs(
+ name='hey',
+ efficiency=None,
+ efficiency_reverse=None,
+ static_loss=[1 for k in range(number_time_intervals)],
+ capacity=tuple(1+o for o in range(number_options)),
+ minimum_cost=tuple(1+o for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ validate=True)
+ except TypeError:
+ error_triggered = True
+ assert error_triggered
+
+ # ValueError(
+ # 'No static loss values were provided. There should be one'+
+ # ' value per option, scenario and time interval.')
+
+ error_triggered = False
+ try:
+ _ = Arcs(
+ name='hey',
+ efficiency=None,
+ efficiency_reverse=None,
+ static_loss={},
+ capacity=tuple(1+o for o in range(number_options)),
+ minimum_cost=tuple(1+o for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ validate=True)
+ except ValueError:
+ error_triggered = True
+ assert error_triggered
+
+ # *************************************************************************
+ # *************************************************************************
+
+ def test_arc_technologies(self):
+
+ # *********************************************************************
+
+ # create arc technology using instantaneous capacities
+
+ number_scenarios = 2
+ number_options = 4
+ number_time_intervals = 3
+
+ efficiency_dict = {
+ (q,k): 0.85
+ for q in range(number_scenarios)
+ for k in range(number_time_intervals)
+ }
+
+ for capacity_is_instantaneous in (True, False):
+
+ arc_tech = Arcs(
+ name='any',
+ efficiency=efficiency_dict,
+ efficiency_reverse=None,
+ static_loss=None,
+ capacity=tuple(1+o for o in range(number_options)),
+ minimum_cost=tuple(1+o for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ validate=True)
+
+ assert arc_tech.has_proportional_losses()
+
+ assert not arc_tech.has_static_losses()
+
+ assert not arc_tech.is_infinite_capacity()
+
+ assert not arc_tech.has_been_selected()
+
+ assert arc_tech.is_isotropic(reverse_none_means_isotropic=True)
+
+ assert not arc_tech.is_isotropic(reverse_none_means_isotropic=False)
+
+ assert arc_tech.has_constant_efficiency()
+
+ # create arc technology with only one option
+
+ arc_tech = Arcs(
+ name='any',
+ efficiency=efficiency_dict,
+ efficiency_reverse=None,
+ static_loss=None,
+ capacity=(1,),
+ minimum_cost=(1,),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ validate=True)
+
+ assert arc_tech.has_proportional_losses()
+
+ assert not arc_tech.has_static_losses()
+
+ assert not arc_tech.is_infinite_capacity()
+
+ assert not arc_tech.has_been_selected()
+
+ assert arc_tech.is_isotropic(reverse_none_means_isotropic=True)
+
+ assert not arc_tech.is_isotropic(reverse_none_means_isotropic=False)
+
+ assert arc_tech.has_constant_efficiency()
+
+ # create arc technology for one time interval
+
+ arc_tech = Arcs(
+ name='any',
+ efficiency={(0,0): 0.95},
+ efficiency_reverse=None,
+ static_loss=None,
+ capacity=tuple(1+o for o in range(number_options)),
+ minimum_cost=tuple(1+o for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ validate=True)
+
+ assert arc_tech.has_proportional_losses()
+
+ assert not arc_tech.has_static_losses()
+
+ assert not arc_tech.is_infinite_capacity()
+
+ assert not arc_tech.has_been_selected()
+
+ assert arc_tech.is_isotropic(reverse_none_means_isotropic=True)
+
+ assert not arc_tech.is_isotropic(reverse_none_means_isotropic=False)
+
+ assert arc_tech.has_constant_efficiency()
+
+ # create arc technology for one time interval and isotropic
+
+ arc_tech = Arcs(
+ name='any',
+ efficiency={(0,0): 0.95},
+ efficiency_reverse={(0,0): 0.95},
+ static_loss=None,
+ capacity=tuple(1+o for o in range(number_options)),
+ minimum_cost=tuple(1+o for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ validate=True)
+
+ assert arc_tech.has_proportional_losses()
+
+ assert not arc_tech.has_static_losses()
+
+ assert not arc_tech.is_infinite_capacity()
+
+ assert not arc_tech.has_been_selected()
+
+ assert arc_tech.is_isotropic(reverse_none_means_isotropic=True)
+
+ assert arc_tech.is_isotropic(reverse_none_means_isotropic=False)
+
+ assert arc_tech.has_constant_efficiency()
+
+ # create arc technology for one time interval and anisotropic
+
+ arc_tech = Arcs(
+ name='any',
+ efficiency={(0,0): 0.95},
+ efficiency_reverse={(0,0): 1},
+ static_loss=None,
+ capacity=tuple(1+o for o in range(number_options)),
+ minimum_cost=tuple(1+o for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ validate=True)
+
+ assert arc_tech.has_proportional_losses()
+
+ assert not arc_tech.has_static_losses()
+
+ assert not arc_tech.is_infinite_capacity()
+
+ assert not arc_tech.has_been_selected()
+
+ assert not arc_tech.is_isotropic(reverse_none_means_isotropic=True)
+
+ assert not arc_tech.is_isotropic(reverse_none_means_isotropic=False)
+
+ assert not arc_tech.has_constant_efficiency()
+
+ # create arc technology for one time interval and anisotropic
+
+ arc_tech = Arcs(
+ name='any',
+ efficiency={(0,0): 1},
+ efficiency_reverse={(0,0): 0.95},
+ static_loss=None,
+ capacity=tuple(1+o for o in range(number_options)),
+ minimum_cost=tuple(1+o for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ validate=True)
+
+ assert arc_tech.has_proportional_losses()
+
+ assert not arc_tech.has_static_losses()
+
+ assert not arc_tech.is_infinite_capacity()
+
+ assert not arc_tech.has_been_selected()
+
+ assert not arc_tech.is_isotropic(reverse_none_means_isotropic=True)
+
+ assert not arc_tech.is_isotropic(reverse_none_means_isotropic=False)
+
+ assert not arc_tech.has_constant_efficiency()
+
+ # create arc technology for one time interval and anisotropic
+
+ arc_tech = Arcs(
+ name='any',
+ efficiency={(0,0): 0.95},
+ efficiency_reverse={(0,0): 0.95},
+ static_loss=None,
+ capacity=tuple(1+o for o in range(number_options)),
+ minimum_cost=tuple(1+o for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ validate=True)
+
+ assert arc_tech.has_proportional_losses()
+
+ assert not arc_tech.has_static_losses()
+
+ assert not arc_tech.is_infinite_capacity()
+
+ assert not arc_tech.has_been_selected()
+
+ assert arc_tech.is_isotropic(reverse_none_means_isotropic=True)
+
+ assert arc_tech.is_isotropic(reverse_none_means_isotropic=False)
+
+ assert arc_tech.has_constant_efficiency()
+
+ # *****************************************************************
+ # *****************************************************************
+
+ # trigger errors
+
+ # TypeError('The name attribute is not hashable.')
+
+ error_triggered = False
+ try:
+ _ = Arcs(
+ name=[1,2,3],
+ efficiency=efficiency_dict,
+ efficiency_reverse=None,
+ static_loss=None,
+ capacity=tuple(1+o for o in range(number_options)),
+ minimum_cost=tuple(1+o for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ validate=True)
+ except TypeError:
+ error_triggered = True
+ assert error_triggered
+
+ #TypeError:The efficiency dict keys must be (scenario, interval) tuples
+
+ error_triggered = False
+ try:
+ _ = Arcs(
+ name='hey',
+ efficiency={k:1 for k in range(number_time_intervals)},
+ efficiency_reverse=None,
+ static_loss=None,
+ capacity=tuple(1+o for o in range(number_options)),
+ minimum_cost=tuple(1+o for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ validate=True)
+ except TypeError:
+ error_triggered = True
+ assert error_triggered
+
+ #ValueError( 'The efficiency dict keys must be tuples of size 2.')
+
+ error_triggered = False
+ try:
+ _ = Arcs(
+ name='hey',
+ efficiency={(k,3,4) :1 for k in range(number_time_intervals)},
+ efficiency_reverse=None,
+ static_loss=None,
+ capacity=tuple(1+o for o in range(number_options)),
+ minimum_cost=tuple(1+o for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ validate=True)
+ except ValueError:
+ error_triggered = True
+ assert error_triggered
+
+ # TypeError(The efficiency should be given as a dict or None.')
+
+ error_triggered = False
+ try:
+ _ = Arcs(
+ name='hey',
+ efficiency=[1 for k in range(number_time_intervals)],
+ efficiency_reverse=None,
+ static_loss=None,
+ capacity=tuple(1+o for o in range(number_options)),
+ minimum_cost=tuple(1+o for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ validate=True)
+ except TypeError:
+ error_triggered = True
+ assert error_triggered
+
+ # TypeError('The reverse efficiency has to match the nominal'+
+ # ' one when there are no proportional losses.')
+
+ error_triggered = False
+ try:
+ _ = Arcs(
+ name='hey',
+ efficiency=None,
+ efficiency_reverse={},
+ static_loss=None,
+ capacity=tuple(1+o for o in range(number_options)),
+ minimum_cost=tuple(1+o for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ validate=True)
+ except TypeError:
+ error_triggered = True
+ assert error_triggered
+
+ # TypeError:'The reverse efficiency should be given as a dict or None.'
+
+ error_triggered = False
+ try:
+ _ = Arcs(
+ name='hey',
+ efficiency=efficiency_dict,
+ efficiency_reverse=[1 for k in range(number_time_intervals)],
+ static_loss=None,
+ capacity=tuple(1+o for o in range(number_options)),
+ minimum_cost=tuple(1+o for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ validate=True)
+ except TypeError:
+ error_triggered = True
+ assert error_triggered
+
+ # ValueError(
+ # 'No efficiency values were provided. There should be '+
+ # 'one value per scenario and time interval.')
+
+ error_triggered = False
+ try:
+ _ = Arcs(
+ name='hey',
+ efficiency=efficiency_dict,
+ efficiency_reverse={},
+ static_loss=None,
+ capacity=tuple(1+o for o in range(number_options)),
+ minimum_cost=tuple(1+o for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ validate=True)
+ except ValueError:
+ error_triggered = True
+ assert error_triggered
+
+ # ValueError: The keys for the efficiency dicts do not match.
+
+ error_triggered = False
+ try:
+ _ = Arcs(
+ name='hey',
+ efficiency=efficiency_dict,
+ efficiency_reverse={
+ (key[1],key[0]): value
+ for key, value in efficiency_dict.items()},
+ static_loss=None,
+ capacity=tuple(1+o for o in range(number_options)),
+ minimum_cost=tuple(1+o for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ validate=True)
+ except ValueError:
+ error_triggered = True
+ assert error_triggered
+
+ # TypeError: Efficiency values must be provided as numeric types.
+
+ error_triggered = False
+ try:
+ _ = Arcs(
+ name='hey',
+ efficiency=efficiency_dict,
+ efficiency_reverse={
+ (key[0],key[1]): str(value)
+ for key, value in efficiency_dict.items()},
+ static_loss=None,
+ capacity=tuple(1+o for o in range(number_options)),
+ minimum_cost=tuple(1+o for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ validate=True)
+ except TypeError:
+ error_triggered = True
+ assert error_triggered
+
+ # ValueError('Efficiency values must be positive.')
+
+ error_triggered = False
+ try:
+ _ = Arcs(
+ name='hey',
+ efficiency=efficiency_dict,
+ efficiency_reverse={
+ (key[0],key[1]): -1
+ for key, value in efficiency_dict.items()},
+ static_loss=None,
+ capacity=tuple(1+o for o in range(number_options)),
+ minimum_cost=tuple(1+o for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ validate=True)
+ except ValueError:
+ error_triggered = True
+ assert error_triggered
+
+ #TypeError('The capacity should be given as a list or tuple.')
+
+ error_triggered = False
+ try:
+ _ = Arcs(
+ name='hey',
+ efficiency=efficiency_dict,
+ efficiency_reverse=None,
+ static_loss=None,
+ capacity={o: 1+o for o in range(number_options)},
+ minimum_cost=tuple(1+o for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ validate=True)
+ except TypeError:
+ error_triggered = True
+ assert error_triggered
+
+ # TypeError: The minimum cost values should be given as a list or tuple
+
+ error_triggered = False
+ try:
+ _ = Arcs(
+ name='hey',
+ efficiency=efficiency_dict,
+ efficiency_reverse=None,
+ static_loss=None,
+ capacity=tuple(1+o for o in range(number_options)),
+ minimum_cost={o: 1+o for o in range(number_options)},
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ validate=True
+ )
+ except TypeError:
+ error_triggered = True
+ assert error_triggered
+
+ # TypeError: The specific capacity cost was not given as a numeric type
+
+ error_triggered = False
+ try:
+ _ = Arcs(
+ name='hey',
+ efficiency=efficiency_dict,
+ efficiency_reverse=None,
+ static_loss=None,
+ capacity=tuple(1+o for o in range(number_options)),
+ minimum_cost=tuple(1+o for o in range(number_options)),
+ specific_capacity_cost=[1],
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ validate=True)
+ except TypeError:
+ error_triggered = True
+ assert error_triggered
+
+ # ValueError:The number of capacity and minimum cost entries must match
+
+ error_triggered = False
+ try:
+ _ = Arcs(
+ name='hey',
+ efficiency=efficiency_dict,
+ efficiency_reverse=None,
+ static_loss=None,
+ capacity=tuple(1+o for o in range(number_options+1)),
+ minimum_cost=tuple(1+o for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ validate=True)
+ except ValueError:
+ error_triggered = True
+ assert error_triggered
+
+ # ValueError: No entries for capacity and minimum cost were provided.
+ # At least one option should be provided.
+
+ error_triggered = False
+ try:
+ _ = Arcs(
+ name='hey',
+ efficiency=efficiency_dict,
+ efficiency_reverse=None,
+ static_loss=None,
+ capacity=tuple(),
+ minimum_cost=tuple(),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ validate=True)
+ except ValueError:
+ error_triggered = True
+ assert error_triggered
+
+ # ValueError: No entries for efficiency were provided. There should be
+ # one entry per time interval.
+
+ error_triggered = False
+ try:
+ _ = Arcs(
+ name='hey',
+ efficiency={},
+ efficiency_reverse=None,
+ static_loss=None,
+ capacity=tuple(1+o for o in range(number_options)),
+ minimum_cost=tuple(1+o for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ validate=True)
+ except ValueError:
+ error_triggered = True
+ assert error_triggered
+
+ # ValueError('The number of efficiency values must match the number of
+ # time intervals.')
+
+ arc_tech = Arcs(
+ name='hey',
+ efficiency={
+ (q,k): 0.85
+ for q in range(number_scenarios)
+ for k in range(number_time_intervals+1)
+ },
+ efficiency_reverse=None,
+ static_loss=None,
+ capacity=tuple(1+o for o in range(number_options)),
+ minimum_cost=tuple(1+o for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ validate=True)
+
+ error_triggered = False
+ try:
+ arc_tech.validate_sizes(number_options=number_options,
+ number_scenarios=number_scenarios,
+ number_intervals=[
+ number_time_intervals
+ for _ in range(number_scenarios)])
+ except ValueError:
+ error_triggered = True
+ assert error_triggered
+
+ # ValueError('The number of efficiency values must match the number of
+ # time intervals.')
+
+ error_triggered = False
+ try:
+ arc_tech = Arcs(
+ name='hey',
+ efficiency={
+ (q,k): 0.85
+ for q in range(number_scenarios)
+ for k in range(number_time_intervals)
+ },
+ efficiency_reverse={
+ (q,k): 0.85
+ for q in range(number_scenarios)
+ for k in range(number_time_intervals-1)
+ },
+ static_loss=None,
+ capacity=tuple(1+o for o in range(number_options)),
+ minimum_cost=tuple(1+o for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ validate=True)
+ arc_tech.validate_sizes(number_options=number_options,
+ number_scenarios=number_scenarios,
+ number_intervals=[
+ number_time_intervals
+ for _ in range(number_scenarios)])
+ except ValueError:
+ error_triggered = True
+ assert error_triggered
+
+ # ValueError('The number of capacity values must match the number of
+ # options.')
+
+ arc_tech = Arcs(
+ name='hey',
+ efficiency=efficiency_dict,
+ efficiency_reverse=None,
+ static_loss=None,
+ capacity=tuple(1+o for o in range(number_options+1)),
+ minimum_cost=tuple(1+o for o in range(number_options+1)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ validate=True
+ )
+
+ error_triggered = False
+ try:
+ arc_tech.validate_sizes(number_options=number_options,
+ number_scenarios=number_scenarios,
+ number_intervals=[
+ number_time_intervals
+ for _ in range(number_scenarios)])
+ except ValueError:
+ error_triggered = True
+ assert error_triggered
+
+ # ValueError: The minimum cost values are inconsistent with the number
+ # of options.
+
+ arc_tech = Arcs(
+ name='hey',
+ efficiency=efficiency_dict,
+ efficiency_reverse=None,
+ static_loss=None,
+ capacity=tuple(1+o for o in range(number_options+1)),
+ minimum_cost=tuple(1+o for o in range(number_options+1)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ validate=True
+ )
+
+ error_triggered = False
+ try:
+ arc_tech.validate_sizes(number_options=number_options,
+ number_scenarios=number_scenarios,
+ number_intervals=[
+ number_time_intervals
+ for _ in range(number_scenarios)])
+ except ValueError:
+ error_triggered = True
+ assert error_triggered
+
+ # TypeError('Efficiency values must be provided as numeric types.')
+
+ error_triggered = False
+ try:
+ _ = Arcs(
+ name='hey',
+ efficiency={key: str(value)
+ for key, value in efficiency_dict.items()},
+ efficiency_reverse=None,
+ static_loss=None,
+ capacity=tuple(1+o for o in range(number_options)),
+ minimum_cost=tuple(1+o for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ validate=True)
+ except TypeError:
+ error_triggered = True
+ assert error_triggered
+
+ # ValueError('Efficiency values must be positive.')
+
+ error_triggered = False
+ try:
+ _ = Arcs(
+ name='hey',
+ efficiency={key: -value*random.randint(0,1)
+ for key, value in efficiency_dict.items()},
+ efficiency_reverse=None,
+ static_loss=None,
+ capacity=tuple(1+o for o in range(number_options)),
+ minimum_cost=tuple(1+o for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ validate=True)
+ except ValueError:
+ error_triggered = True
+ assert error_triggered
+
+ # TypeError('Capacity values must be provided as numeric types.')
+
+ error_triggered = False
+ try:
+ _ = Arcs(
+ name='hey',
+ efficiency=efficiency_dict,
+ efficiency_reverse=None,
+ static_loss=None,
+ capacity=tuple(str(1+o) for o in range(number_options)),
+ minimum_cost=tuple(1+o for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ validate=True)
+ except TypeError:
+ error_triggered = True
+ assert error_triggered
+
+ # ValueError('Capacity values must be positive.')
+
+ error_triggered = False
+ try:
+ _ = Arcs(
+ name='hey',
+ efficiency=efficiency_dict,
+ efficiency_reverse=None,
+ static_loss=None,
+ capacity=tuple(-random.randint(0,1)
+ for o in range(number_options)),
+ minimum_cost=tuple(1+o for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ validate=True)
+ except ValueError:
+ error_triggered = True
+ assert error_triggered
+
+ # TypeError('Minimum cost values must be provided as numeric types.')
+
+ error_triggered = False
+ try:
+ _ = Arcs(
+ name='hey',
+ efficiency=efficiency_dict,
+ efficiency_reverse=None,
+ static_loss=None,
+ capacity=tuple(1+o for o in range(number_options)),
+ minimum_cost=tuple(str(1+o) for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ validate=True)
+ except TypeError:
+ error_triggered = True
+ assert error_triggered
+
+ # ValueError('Minimum cost values must be positive or zero.')
+
+ error_triggered = False
+ try:
+ _ = Arcs(
+ name='hey',
+ efficiency=efficiency_dict,
+ efficiency_reverse=None,
+ static_loss=None,
+ capacity=tuple(1+o
+ for o in range(number_options)),
+ minimum_cost=tuple(-1
+ for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=capacity_is_instantaneous,
+ validate=True)
+ except ValueError:
+ error_triggered = True
+ assert error_triggered
+
+ # TypeError('The information about capacities being instantaneous or not
+ # should be given as a boolean variable.')
+
+ error_triggered = False
+ try:
+ _ = Arcs(
+ name='hey',
+ efficiency=efficiency_dict,
+ efficiency_reverse=None,
+ static_loss=None,
+ capacity=tuple(1+o for o in range(number_options)),
+ minimum_cost=tuple(1+o for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=1,
+ validate=True)
+ except TypeError:
+ error_triggered = True
+ assert error_triggered
+
+ # *********************************************************************
+ # *********************************************************************
+
+ # Network
+
+ arc_tech_AB = Arcs(
+ name='AB',
+ efficiency=efficiency_dict,
+ efficiency_reverse=None,
+ static_loss=None,
+ capacity=tuple(1+o for o in range(number_options)),
+ minimum_cost=tuple(1+o for o in range(number_options)),
+ specific_capacity_cost=1,
+ capacity_is_instantaneous=False,
+ validate=True)
+
+ arc_tech_AB.options_selected[0] = True
+
+ assert arc_tech_AB.number_options() == number_options
+
+ net = Network()
+
+ # add undirected arc
+
+ net.add_undirected_arc(
+ node_key_a='A',
+ node_key_b='B',
+ arcs=arc_tech_AB)
+
+ # add directed arc
+
+ net.add_directed_arc(
+ node_key_a='A',
+ node_key_b='B',
+ arcs=arc_tech_AB)
+
+ # add infinite capacity arc
+
+ net.add_infinite_capacity_arc(
+ node_key_a='C',
+ node_key_b='D',
+ efficiency={
+ (i, j): 1
+ for i in range(3)
+ for j in range(4)},
+ static_loss=None)
+
+ # add pre-existing directed arc
+
+ net.add_preexisting_directed_arc(
+ node_key_a='E',
+ node_key_b='F',
+ efficiency=efficiency_dict,
+ static_loss=None,
+ capacity=3,
+ capacity_is_instantaneous=True)
+
+ # add pre-existing undirected arc
+
+ net.add_preexisting_undirected_arc(
+ node_key_a='A',
+ node_key_b='C',
+ efficiency=efficiency_dict,
+ efficiency_reverse=efficiency_dict,
+ static_loss=None,
+ capacity=3,
+ capacity_is_instantaneous=True)
+
+ net.modify_network_arc(
+ node_key_a='A',
+ node_key_b='C',
+ arc_key_ab='AC',
+ data_dict={net.KEY_ARC_TECH: arc_tech_AB, net.KEY_ARC_UND: False})
+
+ # *********************************************************************
+ # *********************************************************************
+
+ # add import node
+
+ imp_resource_price = ResourcePrice(
+ prices=[random.random()
+ for k in range(number_time_intervals)],
+ volumes=[ *[random.random() for k in range(number_time_intervals-1)], None]
+ )
+
+ net.add_import_node(node_key='G', prices={(0,0,0): imp_resource_price})
+
+ # add export node
+
+ exp_resource_price = ResourcePrice(
+ prices=[random.random()
+ for k in range(number_time_intervals)],
+ volumes=[ *[random.random() for k in range(number_time_intervals-1)], None]
+ )
+
+ net.add_export_node(node_key='H', prices={(0,0,0): exp_resource_price})
+
+ net.add_waypoint_node(node_key='Z')
+
+ base_flow = {
+ (i,j): random.random()
+ for i in range(3)
+ for j in range(4)
+ }
+
+ net.add_source_sink_node(node_key='Y', base_flow=base_flow)
+
+ base_flow[(2,3)] = random.random()
+
+ net.modify_network_node(
+ node_key='Y',
+ node_data={net.KEY_NODE_BASE_FLOW: base_flow}
+ )
+
+ net.identify_node_types()
+
+ assert 'Z' in net.waypoint_nodes
+
+ assert 'G' in net.import_nodes
+
+ assert 'H' in net.export_nodes
+
+ assert 'Y' in net.source_sink_nodes
+
+ # *************************************************************************
+ # *************************************************************************
+
+ def test_arcs_without_losses(self):
+
+ # test arc without (static and proportional) losses
+
+ arc_tech = ArcsWithoutLosses(
+ name='AB',
+ capacity=(1,2,3),
+ minimum_cost=(4,5,6),
+ specific_capacity_cost=6,
+ capacity_is_instantaneous=False,
+ validate=True
+ )
+
+ assert not arc_tech.has_proportional_losses()
+
+ assert not arc_tech.has_static_losses()
+
+ assert not arc_tech.is_infinite_capacity()
+
+ assert arc_tech.has_constant_efficiency()
+
+ # test arc without static losses
+
+ arc_tech = ArcsWithoutStaticLosses(
+ name='AB',
+ efficiency={(0,0):1, (0,1):0.9, (0,2):0.8},
+ efficiency_reverse=None,
+ capacity=(1,2,3),
+ minimum_cost=(4,5,6),
+ specific_capacity_cost=6,
+ capacity_is_instantaneous=False,
+ validate=True
+ )
+
+ assert arc_tech.has_proportional_losses()
+
+ assert not arc_tech.has_static_losses()
+
+ assert not arc_tech.is_infinite_capacity()
+
+ assert not arc_tech.has_constant_efficiency()
+
+ # test arc without proportional losses
+
+ arc_tech = ArcsWithoutProportionalLosses(
+ name='AB',
+ static_loss={(0,0,0):0.1, (0,0,1):0.2, (0,0,2):0.3,
+ (1,0,0):0.15, (1,0,1):0.25, (1,0,2):0.35,
+ (2,0,0):0.16, (2,0,1):0.26, (2,0,2):0.36},
+ capacity=(1,2,3),
+ minimum_cost=(4,5,6),
+ specific_capacity_cost=6,
+ capacity_is_instantaneous=False,
+ validate=True
+ )
+
+ assert not arc_tech.has_proportional_losses()
+
+ assert arc_tech.has_static_losses()
+
+ assert not arc_tech.is_infinite_capacity()
+
+ assert arc_tech.has_constant_efficiency()
+
+ # *************************************************************************
+ # *************************************************************************
+
+ def test_modifying_nodes(self):
+
+ # *********************************************************************
+
+ net = Network()
+
+ number_intervals = 3
+
+ resource_price = ResourcePrice(
+ prices=[random.random() for k in range(number_intervals)],
+ volumes=[
+ *[random.random() for k in range(number_intervals-1)], None
+ ]
+ )
+
+ base_flow = {
+ (0,k): random.random()
+ for k in range(number_intervals)}
+
+ arc_tech = ArcsWithoutLosses(
+ name='hello',
+ capacity=[5],
+ minimum_cost=[3],
+ specific_capacity_cost=3,
+ capacity_is_instantaneous=False
+ )
+
+ # add isolated import node
+
+ net.add_import_node(node_key='I_iso',
+ prices={(0,0,0): resource_price})
+
+ # add import node with outgoing arcs
+
+ net.add_import_node(node_key='I',
+ prices={(0,0,0): resource_price})
+
+ # add isolated export node
+
+ net.add_import_node(node_key='E_iso',
+ prices={(0,0,0): resource_price})
+
+ # add export node with incoming arcs
+
+ net.add_export_node(node_key='E',
+ prices={(0,0,0): resource_price})
+
+ # add isolated normal node
+
+ net.add_source_sink_node(node_key='A_iso',
+ base_flow=base_flow)
+
+ # add normal node with incoming arcs
+
+ net.add_source_sink_node(node_key='A_in',
+ base_flow=base_flow)
+
+ # add normal node with outgoing arcs
+
+ net.add_source_sink_node(node_key='A_out',
+ base_flow=base_flow)
+
+ # add normal node with incoming and outgoing arcs
+
+ net.add_source_sink_node(node_key='A',
+ base_flow=base_flow)
+
+ # *********************************************************************
+
+ # arcs
+
+ net.add_directed_arc(node_key_a='I',
+ node_key_b='A_in',
+ arcs=arc_tech)
+
+ net.add_directed_arc(node_key_a='I',
+ node_key_b='A',
+ arcs=arc_tech)
+
+ net.add_directed_arc(node_key_a='A_out',
+ node_key_b='E',
+ arcs=arc_tech)
+
+ net.add_directed_arc(node_key_a='A',
+ node_key_b='E',
+ arcs=arc_tech)
+
+ # *********************************************************************
+
+ # change I_iso to regular: okay
+
+ net.modify_network_node(
+ node_key='I_iso',
+ node_data={net.KEY_NODE_TYPE: net.KEY_NODE_TYPE_SOURCE_SINK,
+ net.KEY_NODE_BASE_FLOW: base_flow})
+
+ # reverse: okay
+
+ net.modify_network_node(
+ node_key='I_iso',
+ node_data={net.KEY_NODE_TYPE: net.KEY_NODE_TYPE_IMP,
+ net.KEY_NODE_PRICES: resource_price})
+
+ # change I_iso to export: okay
+
+ net.modify_network_node(
+ node_key='I_iso',
+ node_data={net.KEY_NODE_TYPE: net.KEY_NODE_TYPE_EXP,
+ net.KEY_NODE_PRICES: resource_price})
+
+ # reverse: okay
+
+ net.modify_network_node(
+ node_key='I_iso',
+ node_data={net.KEY_NODE_TYPE: net.KEY_NODE_TYPE_IMP,
+ net.KEY_NODE_PRICES: resource_price})
+
+ # change I_iso to waypoint: okay
+
+ net.modify_network_node(
+ node_key='I_iso',
+ node_data={net.KEY_NODE_TYPE: net.KEY_NODE_TYPE_WAY})
+
+ # reverse: okay
+
+ net.modify_network_node(
+ node_key='I_iso',
+ node_data={net.KEY_NODE_TYPE: net.KEY_NODE_TYPE_IMP,
+ net.KEY_NODE_PRICES: resource_price})
+
+ # *********************************************************************
+
+ # change E_iso to regular: okay
+
+ net.modify_network_node(
+ node_key='E_iso',
+ node_data={net.KEY_NODE_TYPE: net.KEY_NODE_TYPE_SOURCE_SINK,
+ net.KEY_NODE_BASE_FLOW: base_flow})
+
+ # reverse: okay
+
+ net.modify_network_node(
+ node_key='E_iso',
+ node_data={net.KEY_NODE_TYPE: net.KEY_NODE_TYPE_EXP,
+ net.KEY_NODE_PRICES: resource_price})
+
+ # change E_iso to import: okay
+
+ net.modify_network_node(
+ node_key='E_iso',
+ node_data={net.KEY_NODE_TYPE: net.KEY_NODE_TYPE_IMP,
+ net.KEY_NODE_PRICES: resource_price})
+
+ # reverse: okay
+
+ net.modify_network_node(
+ node_key='E_iso',
+ node_data={net.KEY_NODE_TYPE: net.KEY_NODE_TYPE_EXP,
+ net.KEY_NODE_PRICES: resource_price})
+
+ # change E_iso to waypoint: okay
+
+ net.modify_network_node(
+ node_key='E_iso',
+ node_data={net.KEY_NODE_TYPE: net.KEY_NODE_TYPE_WAY})
+
+ # reverse: okay
+
+ net.modify_network_node(
+ node_key='E_iso',
+ node_data={net.KEY_NODE_TYPE: net.KEY_NODE_TYPE_EXP,
+ net.KEY_NODE_PRICES: resource_price})
+
+ # *********************************************************************
+
+ # change A_iso to export: okay
+
+ net.modify_network_node(
+ node_key='A_iso',
+ node_data={net.KEY_NODE_TYPE: net.KEY_NODE_TYPE_EXP,
+ net.KEY_NODE_PRICES: resource_price})
+
+ # reverse: okay
+
+ net.modify_network_node(
+ node_key='A_iso',
+ node_data={net.KEY_NODE_TYPE: net.KEY_NODE_TYPE_SOURCE_SINK,
+ net.KEY_NODE_BASE_FLOW: base_flow})
+
+ # change A_iso to import: okay
+
+ net.modify_network_node(
+ node_key='A_iso',
+ node_data={net.KEY_NODE_TYPE: net.KEY_NODE_TYPE_IMP,
+ net.KEY_NODE_PRICES: resource_price})
+
+ # reverse: okay
+
+ net.modify_network_node(
+ node_key='A_iso',
+ node_data={net.KEY_NODE_TYPE: net.KEY_NODE_TYPE_SOURCE_SINK,
+ net.KEY_NODE_BASE_FLOW: base_flow})
+
+ # change A_iso to waypoint: okay
+
+ net.modify_network_node(
+ node_key='A_iso',
+ node_data={net.KEY_NODE_TYPE: net.KEY_NODE_TYPE_WAY})
+
+ # reverse: okay
+
+ net.modify_network_node(
+ node_key='A_iso',
+ node_data={net.KEY_NODE_TYPE: net.KEY_NODE_TYPE_SOURCE_SINK,
+ net.KEY_NODE_BASE_FLOW: base_flow})
+
+ # *********************************************************************
+
+ # change I to regular: okay
+
+ net.modify_network_node(
+ node_key='I',
+ node_data={net.KEY_NODE_TYPE: net.KEY_NODE_TYPE_SOURCE_SINK,
+ net.KEY_NODE_BASE_FLOW: base_flow})
+
+ # reverse: okay
+
+ net.modify_network_node(
+ node_key='I',
+ node_data={net.KEY_NODE_TYPE: net.KEY_NODE_TYPE_IMP,
+ net.KEY_NODE_PRICES: resource_price})
+
+ # change I to waypoint: okay
+
+ net.modify_network_node(
+ node_key='I',
+ node_data={net.KEY_NODE_TYPE: net.KEY_NODE_TYPE_WAY})
+
+ # reverse: okay
+
+ net.modify_network_node(
+ node_key='I',
+ node_data={net.KEY_NODE_TYPE: net.KEY_NODE_TYPE_IMP,
+ net.KEY_NODE_PRICES: resource_price})
+
+ # *********************************************************************
+
+ # change E to regular: okay
+
+ net.modify_network_node(
+ node_key='E',
+ node_data={net.KEY_NODE_TYPE: net.KEY_NODE_TYPE_SOURCE_SINK,
+ net.KEY_NODE_BASE_FLOW: base_flow})
+
+ # reverse: okay
+
+ net.modify_network_node(
+ node_key='E',
+ node_data={net.KEY_NODE_TYPE: net.KEY_NODE_TYPE_EXP,
+ net.KEY_NODE_PRICES: resource_price})
+
+ # change E to waypoint: okay
+
+ net.modify_network_node(
+ node_key='E',
+ node_data={net.KEY_NODE_TYPE: net.KEY_NODE_TYPE_WAY})
+
+ # reverse: okay
+
+ net.modify_network_node(
+ node_key='E',
+ node_data={net.KEY_NODE_TYPE: net.KEY_NODE_TYPE_EXP,
+ net.KEY_NODE_PRICES: resource_price})
+
+ # *********************************************************************
+
+ # change A_in to export: okay
+
+ net.modify_network_node(
+ node_key='A_in',
+ node_data={net.KEY_NODE_TYPE: net.KEY_NODE_TYPE_EXP,
+ net.KEY_NODE_PRICES: resource_price})
+
+ # reverse: okay
+
+ net.modify_network_node(
+ node_key='A_in',
+ node_data={net.KEY_NODE_TYPE: net.KEY_NODE_TYPE_SOURCE_SINK,
+ net.KEY_NODE_BASE_FLOW: base_flow})
+
+ # change A_in to waypoint: okay
+
+ net.modify_network_node(
+ node_key='A_in',
+ node_data={net.KEY_NODE_TYPE: net.KEY_NODE_TYPE_WAY})
+
+ # reverse: okay
+
+ net.modify_network_node(
+ node_key='A_in',
+ node_data={net.KEY_NODE_TYPE: net.KEY_NODE_TYPE_SOURCE_SINK,
+ net.KEY_NODE_BASE_FLOW: base_flow})
+
+ # *********************************************************************
+
+ # change A_out to import: okay
+
+ net.modify_network_node(
+ node_key='A_out',
+ node_data={net.KEY_NODE_TYPE: net.KEY_NODE_TYPE_IMP,
+ net.KEY_NODE_PRICES: resource_price})
+
+ # reverse: okay
+
+ net.modify_network_node(
+ node_key='A_out',
+ node_data={net.KEY_NODE_TYPE: net.KEY_NODE_TYPE_SOURCE_SINK,
+ net.KEY_NODE_BASE_FLOW: base_flow})
+
+ # change A_out to waypoint: okay
+
+ net.modify_network_node(
+ node_key='A_out',
+ node_data={net.KEY_NODE_TYPE: net.KEY_NODE_TYPE_WAY})
+
+ # reverse: okay
+
+ net.modify_network_node(
+ node_key='A_out',
+ node_data={net.KEY_NODE_TYPE: net.KEY_NODE_TYPE_SOURCE_SINK,
+ net.KEY_NODE_BASE_FLOW: base_flow})
+
+ # *********************************************************************
+
+ # change I to export: fail
+
+ error_triggered = False
+ try:
+ net.modify_network_node(
+ node_key='I',
+ node_data={net.KEY_NODE_TYPE: net.KEY_NODE_TYPE_EXP,
+ net.KEY_NODE_PRICES: resource_price}
+ )
+ except ValueError:
+ error_triggered = True
+ assert error_triggered
+
+ # change E to import: fail
+
+ error_triggered = False
+ try:
+ net.modify_network_node(
+ node_key='E',
+ node_data={net.KEY_NODE_TYPE: net.KEY_NODE_TYPE_IMP,
+ net.KEY_NODE_PRICES: resource_price}
+ )
+ except ValueError:
+ error_triggered = True
+ assert error_triggered
+
+ # change A_out to export: fail
+
+ error_triggered = False
+ try:
+ net.modify_network_node(
+ node_key='A_out',
+ node_data={net.KEY_NODE_TYPE: net.KEY_NODE_TYPE_EXP,
+ net.KEY_NODE_PRICES: resource_price}
+ )
+ except ValueError:
+ error_triggered = True
+ assert error_triggered
+
+ # change A_in to import: fail
+
+ error_triggered = False
+ try:
+ net.modify_network_node(
+ node_key='A_in',
+ node_data={net.KEY_NODE_TYPE: net.KEY_NODE_TYPE_IMP,
+ net.KEY_NODE_PRICES: resource_price}
+ )
+ except ValueError:
+ error_triggered = True
+ assert error_triggered
+
+ # change A to export: fail
+
+ error_triggered = False
+ try:
+ net.modify_network_node(
+ node_key='A',
+ node_data={net.KEY_NODE_TYPE: net.KEY_NODE_TYPE_EXP,
+ net.KEY_NODE_PRICES: resource_price}
+ )
+ except ValueError:
+ error_triggered = True
+ assert error_triggered
+
+ # change A to import: fail
+
+ error_triggered = False
+ try:
+ net.modify_network_node(
+ node_key='A',
+ node_data={net.KEY_NODE_TYPE: net.KEY_NODE_TYPE_IMP,
+ net.KEY_NODE_PRICES: resource_price}
+ )
+ except ValueError:
+ error_triggered = True
+ assert error_triggered
+
+ # *********************************************************************
+
+ # try to modify a non-existent node
+
+ error_triggered = False
+ try:
+ net.modify_network_node(
+ node_key='ABCD',
+ node_data={net.KEY_NODE_TYPE: net.KEY_NODE_TYPE_WAY}
+ )
+ except ValueError:
+ error_triggered = True
+ assert error_triggered
+
+ # *********************************************************************
+
+ # *************************************************************************
+ # *************************************************************************
+
+ def test_network_disallowed_cases(self):
+
+ # *********************************************************************
+
+ net = Network()
+
+ number_intervals = 3
+
+ resource_price = ResourcePrice(
+ prices=[random.random() for k in range(number_intervals)],
+ volumes=[
+ *[random.random() for k in range(number_intervals-1)], None
+ ]
+ )
+
+ base_flow = {
+ (0,k): random.random()
+ for k in range(number_intervals)}
+
+ lossless_arcs = ArcsWithoutLosses(
+ name='hello',
+ capacity=[5],
+ minimum_cost=[3],
+ specific_capacity_cost=3,
+ capacity_is_instantaneous=False
+ )
+
+ lossy_arcs = ArcsWithoutProportionalLosses(
+ name='hello back',
+ static_loss={
+ (0,0,k): random.random()
+ for k in range(number_intervals)
+ },
+ capacity=(1,),
+ minimum_cost=(5,),
+ specific_capacity_cost=0,
+ capacity_is_instantaneous=False
+ )
+
+ # add import node I
+
+ net.add_import_node(node_key='I',
+ prices={(0,0,0): resource_price})
+
+ # add export node E
+
+ net.add_export_node(node_key='E',
+ prices={(0,0,0): resource_price})
+
+ # add regular node A
+
+ net.add_source_sink_node(node_key='A',
+ base_flow=base_flow)
+
+ # add regular node B
+
+ net.add_source_sink_node(node_key='B',
+ base_flow=base_flow)
+
+ # add a valid import-export arc
+
+ net.add_directed_arc(node_key_a='I',
+ node_key_b='E',
+ arcs=lossless_arcs)
+
+ # identify the nodes and validate
+
+ net.identify_node_types()
+
+ # *********************************************************************
+ # *********************************************************************
+
+ # trigger errors using pre-identified nodes
+
+ # directed arcs cannot start in an export node: E -> B
+
+ error_triggered = False
+ try:
+ net.add_directed_arc(node_key_a='E',
+ node_key_b='B',
+ arcs=lossless_arcs)
+ except ValueError:
+ error_triggered = True
+ assert error_triggered
+
+ # directed arcs cannot end on an import node: A -> I
+
+ error_triggered = False
+ try:
+ net.add_directed_arc(node_key_a='A',
+ node_key_b='I',
+ arcs=lossless_arcs)
+ except ValueError:
+ error_triggered = True
+ assert error_triggered
+
+ # import-export nodes cannot have static losses
+
+ error_triggered = False
+ try:
+ net.add_directed_arc(node_key_a='I',
+ node_key_b='E',
+ arcs=lossy_arcs)
+ except ValueError:
+ error_triggered = True
+ assert error_triggered
+
+ # undirected arcs cannot involve import nor export nodes
+
+ error_triggered = False
+ try:
+ net.add_undirected_arc(node_key_a='I',
+ node_key_b='A',
+ arcs=lossless_arcs)
+ except ValueError:
+ error_triggered = True
+ assert error_triggered
+
+ # undirected arcs cannot involve import nor export nodes
+
+ error_triggered = False
+ try:
+ net.add_undirected_arc(node_key_a='B',
+ node_key_b='E',
+ arcs=lossless_arcs)
+ except ValueError:
+ error_triggered = True
+ assert error_triggered
+
+ # *********************************************************************
+ # *********************************************************************
+
+ # trigger errors using non-identified nodes
+
+ # *********************************************************************
+
+ # create a new export node
+
+ net.add_export_node(node_key='E1',
+ prices={(0,0,0): resource_price})
+
+ # create an arc starting in that export node
+
+ error_triggered = False
+ try:
+ net.add_directed_arc(node_key_a='E1',
+ node_key_b='B',
+ arcs=lossless_arcs)
+ net.identify_node_types()
+ except ValueError:
+ error_triggered = True
+ assert error_triggered
+
+ # remove the troublesome arc
+
+ net.remove_edge(u='E1', v='B')
+
+ # *********************************************************************
+
+ # create a new import node
+
+ net.add_import_node(node_key='I1',
+ prices={(0,0,0): resource_price})
+
+ # create an arc ending in that import node
+
+ error_triggered = False
+ try:
+ net.add_directed_arc(node_key_a='A',
+ node_key_b='I1',
+ arcs=lossless_arcs)
+ net.identify_node_types()
+ except ValueError:
+ error_triggered = True
+ assert error_triggered
+
+ # remove the troublesome arc
+
+ net.remove_edge(u='A', v='I1')
+
+ # *********************************************************************
+
+ # check non-existent arc
+
+ net.arc_is_undirected(('X','Y', 1))
+
+ # *************************************************************************
+ # *************************************************************************
+
+ def test_pseudo_unique_key_generation(self):
+
+ # create network
+
+ network = Network()
+
+ # add node A
+
+ network.add_waypoint_node(node_key='A')
+
+ # add node B
+
+ network.add_waypoint_node(node_key='B')
+
+ # identify nodes
+
+ network.identify_node_types()
+
+ # add arcs
+
+ key_list = ['3e225573-4e78-48c8-bb08-efbeeb795c22',
+ 'f6d30428-15d1-41e9-a952-0742eaaa5a31',
+ '8c29b906-2518-41c5-ada8-07b83508b5b8',
+ 'f9a72a39-1422-4a02-af97-906ce79c32a3',
+ 'b6941a48-10cc-465d-bf53-178bd2939bd1']
+
+ for key in key_list:
+
+ network.add_edge(
+ u_for_edge='A',
+ v_for_edge='B',
+ key=key,
+ **{network.KEY_ARC_UND: False,
+ network.KEY_ARC_TECH: None}
+ )
+
+ # use a seed number to trigger more iterations
+
+ import uuid
+ rand = random.Random()
+ rand.seed(360)
+ uuid.uuid4 = lambda: uuid.UUID(int=rand.getrandbits(128), version=4)
+
+ error_triggered = False
+ try:
+ _ = network.get_pseudo_unique_arc_key(
+ node_key_start='A',
+ node_key_end='B',
+ max_iterations=len(key_list)-1)
+ except Exception:
+ error_triggered = True
+ assert error_triggered
+
+# *****************************************************************************
+# *****************************************************************************
\ No newline at end of file
diff --git a/tests/test_esipp_problem.py b/tests/test_esipp_problem.py
index 26d15f9..862e5cf 100644
--- a/tests/test_esipp_problem.py
+++ b/tests/test_esipp_problem.py
@@ -569,7 +569,35 @@ class TestESIPPProblem:
# the objective function should be -9.7
assert math.isclose(pyo.value(ipp.instance.obj_f), -9.7, abs_tol=1e-3)
- # TODO: compare model pretty print with expected one
+ # TODO: create method to automate getting data from the command line
+ import io
+ import sys
+ from contextlib import redirect_stdout
+ # print('wow wow wow')
+ # ipp.instance.constr_imp_flow_cost.pprint()
+ expected_string = """constr_imp_flow_cost : Size=4, Index=constr_imp_flow_cost_index, Active=True\n Key : Lower : Body : Upper : Active\n ('mynet', 'thatimpnode', 'peak', 0, 0) : 0.0 : 0*var_if_glqpks[mynet,thatimpnode,peak,0,0,0] - var_ifc_glqpk[mynet,thatimpnode,peak,0,0] : 0.0 : True\n ('mynet', 'thatimpnode', 'peak', 1, 0) : 0.0 : 0*var_if_glqpks[mynet,thatimpnode,peak,1,0,0] - var_ifc_glqpk[mynet,thatimpnode,peak,1,0] : 0.0 : True\n ('mynet', 'thatimpnode', 'total', 0, 0) : 0.0 : var_if_glqpks[mynet,thatimpnode,total,0,0,0] - var_ifc_glqpk[mynet,thatimpnode,total,0,0] : 0.0 : True\n ('mynet', 'thatimpnode', 'total', 1, 0) : 0.0 : var_if_glqpks[mynet,thatimpnode,total,1,0,0] - var_ifc_glqpk[mynet,thatimpnode,total,1,0] : 0.0 : True\n"""
+
+ cmd_output = io.StringIO()
+ sys.stdout = cmd_output
+ ipp.instance.constr_imp_flow_cost.pprint()
+ sys.stdout = sys.__stdout__
+ assert cmd_output.getvalue() == expected_string
+
+ expected_string = """constr_exp_flow_revenue : Size=0, Index=constr_exp_flow_revenue_index, Active=True\n Key : Lower : Body : Upper : Active\n"""
+ f = io.StringIO()
+ with redirect_stdout(f):
+ ipp.instance.constr_exp_flow_revenue.pprint()
+ assert f.getvalue() == expected_string
+
+ # try the whole model
+ # print('wow wow wow')
+ # ipp.instance.pprint()
+ # expected_string = """constr_imp_flow_cost : Size=4, Index=constr_imp_flow_cost_index, Active=True\n Key : Lower : Body : Upper : Active\n ('mynet', 'thatimpnode', 'peak', 0, 0) : 0.0 : 0*var_if_glqpks[mynet,thatimpnode,peak,0,0,0] - var_ifc_glqpk[mynet,thatimpnode,peak,0,0] : 0.0 : True\n ('mynet', 'thatimpnode', 'peak', 1, 0) : 0.0 : 0*var_if_glqpks[mynet,thatimpnode,peak,1,0,0] - var_ifc_glqpk[mynet,thatimpnode,peak,1,0] : 0.0 : True\n ('mynet', 'thatimpnode', 'total', 0, 0) : 0.0 : var_if_glqpks[mynet,thatimpnode,total,0,0,0] - var_ifc_glqpk[mynet,thatimpnode,total,0,0] : 0.0 : True\n ('mynet', 'thatimpnode', 'total', 1, 0) : 0.0 : var_if_glqpks[mynet,thatimpnode,total,1,0,0] - var_ifc_glqpk[mynet,thatimpnode,total,1,0] : 0.0 : True\n"""
+ # cmd_output = io.StringIO()
+ # sys.stdout = cmd_output
+ # ipp.instance.pprint()
+ # sys.stdout = sys.__stdout__
+ # assert cmd_output.getvalue() == expected_string
# from contextlib import redirect_stdout
# import io
--
GitLab