From 62e6ad6d7f0182b9ebbcb015328c90f75f447e21 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Pedro=20L=2E=20Magalh=C3=A3es?= <pmlpm@posteo.de>
Date: Sun, 12 May 2024 00:35:16 +0200
Subject: [PATCH] Added method to check for overlapping segments.
---
src/topupopt/data/gis/utils.py | 83 ++++++++++++++++--
tests/test_esipp_problem.py | 151 ++++++++++++++++++++++++++++++++-
tests/test_gis_utils.py | 24 +++++-
3 files changed, 248 insertions(+), 10 deletions(-)
diff --git a/src/topupopt/data/gis/utils.py b/src/topupopt/data/gis/utils.py
index 0c21992..4197447 100644
--- a/src/topupopt/data/gis/utils.py
+++ b/src/topupopt/data/gis/utils.py
@@ -13,6 +13,7 @@ from numpy import float64, int64
from geopandas import GeoDataFrame, read_file
from shapely.geometry import Point, LineString
import contextily as cx
+from shapely import intersects
# local, internal
@@ -1232,15 +1233,87 @@ def convert_edge_path(
# *****************************************************************************
# *****************************************************************************
-def create_edge_geometry(network: MultiDiGraph, edge_key) -> LineString:
+def create_edge_geometry(
+ network: MultiDiGraph,
+ edge_key,
+ x_key = osm.KEY_OSMNX_X,
+ y_key = osm.KEY_OSMNX_Y) -> LineString:
"Returns a newly-created geometry for a given edge."
return LineString(
- [(network.nodes[edge_key[0]][osm.KEY_OSMNX_X],
- network.nodes[edge_key[0]][osm.KEY_OSMNX_Y]),
- (network.nodes[edge_key[1]][osm.KEY_OSMNX_X],
- network.nodes[edge_key[1]][osm.KEY_OSMNX_Y])]
+ [(network.nodes[edge_key[0]][x_key],
+ network.nodes[edge_key[0]][y_key]),
+ (network.nodes[edge_key[1]][x_key],
+ network.nodes[edge_key[1]][y_key])]
)
# *****************************************************************************
# *****************************************************************************
+
+def find_overlapping_edges(
+ network: MultiDiGraph,
+ excluded_edges: list = None
+ ) -> list:
+ """
+ Returns a list of key pairs for edges whose geometries overlap.
+
+ Parameters
+ ----------
+ network : MultiDiGraph
+ The object describing the network.
+ excluded_edges : list, optional
+ A list of edges that should not be considered. The default is None, in
+ which case all edges in the network object will be considered.
+
+ Returns
+ -------
+ list
+ A list containing key pairs for overlapping edges.
+
+ """
+ # check if there are excluded edges
+ if type(excluded_edges) == type(None):
+ excluded_edges = list()
+ # initialise the list of edges to check
+ edges = list(
+ edge_key
+ for edge_key in network.edges(keys=True)
+ if edge_key not in excluded_edges
+ )
+ visited_edges = []
+ out = []
+ # for each edge
+ for edge_key in edges:
+ # remove the current edge so it is not considered again
+ visited_edges.append(edge_key)
+ # for each other edge
+ for other_edge_key in edges:
+ # for each other edge
+ # skip edges having nodes in common
+ # this will also skip identical edges
+ if edge_key[0] in other_edge_key[0:2] or edge_key[1] in other_edge_key[0:2]:
+ # has nodes in common, skip
+ continue
+ # skip edges that have already been considered in the first loop
+ if other_edge_key in visited_edges:
+ # this edge has already been tested against all other edges
+ continue
+ # first edge
+ if osm.KEY_OSMNX_GEOMETRY in network.edges[edge_key]:
+ first_geo = network.edges[edge_key][osm.KEY_OSMNX_GEOMETRY]
+ else:
+ first_geo = create_edge_geometry(network, edge_key)
+ # second edge
+ if osm.KEY_OSMNX_GEOMETRY in network.edges[other_edge_key]:
+ second_geo = network.edges[other_edge_key][osm.KEY_OSMNX_GEOMETRY]
+ else:
+ second_geo = create_edge_geometry(network, other_edge_key)
+ # check if they intersect
+ if intersects(first_geo, second_geo):
+ # they do, add tuple of the edges to the output
+ out.append((edge_key, other_edge_key))
+ # return tuples of overlapping edges
+ return out
+
+# *****************************************************************************
+# *****************************************************************************
diff --git a/tests/test_esipp_problem.py b/tests/test_esipp_problem.py
index eb71285..c69bbc9 100644
--- a/tests/test_esipp_problem.py
+++ b/tests/test_esipp_problem.py
@@ -3639,8 +3639,6 @@ class TestESIPPProblem:
# *************************************************************************
# *************************************************************************
-
- # TODO: trigger error with static losses
def test_direct_imp_exp_network(self):
@@ -8153,7 +8151,7 @@ class TestESIPPProblem:
solver_options={},
perform_analysis=False,
plot_results=False, # True,
- print_solver_output=True,
+ print_solver_output=False,
time_frame=tf,
networks={"mynet": mynet},
static_losses_mode=True, # just to reach a line,
@@ -8196,6 +8194,153 @@ class TestESIPPProblem:
assert math.isclose(pyo.value(ipp.instance.var_capex), 10.0, abs_tol=1e-3)
# the objective function
assert math.isclose(pyo.value(ipp.instance.obj_f), -1.193236715e+01, abs_tol=1e-3)
+
+
+ # *************************************************************************
+ # *************************************************************************
+
+ def test_problem_with_reverse_tree_network(self):
+
+ # assessment
+ q = 0
+ tf = EconomicTimeFrame(
+ discount_rate=3.5/100,
+ reporting_periods={q: (0,)},
+ reporting_period_durations={q: (365 * 24 * 3600,)},
+ time_intervals={q: (0,)},
+ time_interval_durations={q: (1,)},
+ )
+
+ # 2 nodes: one import, one regular
+ mynet = Network()
+
+ # export node
+ node_EXP = "thatexpnode"
+ mynet.add_export_node(
+ node_key=node_EXP,
+ prices={
+ qpk: ResourcePrice(prices=1.0, volumes=None)
+ for qpk in tf.qpk()
+ },
+ )
+
+ # node A
+ node_A = "thatnodea"
+ mynet.add_source_sink_node(
+ node_key=node_A,
+ base_flow={(q, 0): -0.50},
+ )
+ # node B
+ node_B = "thatnodeb"
+ mynet.add_source_sink_node(
+ node_key=node_B,
+ base_flow={(q, 0): -0.25},
+ )
+ # node C
+ node_C = "thatnodec"
+ mynet.add_source_sink_node(
+ node_key=node_C,
+ base_flow={(q, 0): -1.25},
+ )
+
+ list_imp_arcs = [
+ (node_A, node_EXP), # AE
+ (node_B, node_EXP), # BE
+ (node_C, node_EXP), # CE
+ ]
+ for i, node_pair in enumerate(list_imp_arcs):
+
+ # import arcs: AE, BE, CE
+ new_arc = Arcs(
+ name="arc_"+str(node_pair),
+ efficiency=None,
+ efficiency_reverse=None,
+ static_loss=None,
+ capacity=[2],
+ minimum_cost=[6],
+ specific_capacity_cost=i,
+ capacity_is_instantaneous=False,
+ validate=False,
+ )
+ mynet.add_directed_arc(*node_pair, arcs=new_arc)
+
+ # arcs: AB, BA, BC, CB, AC, CA
+
+ list_other_arcs = [
+ (node_A, node_B), # AB
+ (node_B, node_A), # BA
+ (node_B, node_C), # BC
+ (node_C, node_B), # CB
+ (node_A, node_C), # AC
+ (node_C, node_A), # CA
+ ]
+
+ for node_pair in list_other_arcs:
+ # arc
+ new_arc_tech = Arcs(
+ name="any",
+ efficiency=None,
+ efficiency_reverse=None,
+ static_loss=None,
+ capacity=[3],
+ minimum_cost=[2],
+ specific_capacity_cost=0,
+ capacity_is_instantaneous=False,
+ validate=False,
+ )
+ mynet.add_directed_arc(*node_pair, arcs=new_arc_tech)
+
+ # identify node types
+ mynet.identify_node_types()
+
+ # no sos, regular time intervals
+ ipp = self.build_solve_ipp(
+ solver_options={},
+ perform_analysis=False,
+ plot_results=False, # True,
+ print_solver_output=False,
+ time_frame=tf,
+ networks={"mynet": mynet},
+ static_losses_mode=True, # just to reach a line,
+ mandatory_arcs=[],
+ max_number_parallel_arcs={},
+ simplify_problem=True,
+ )
+
+ assert ipp.has_peak_total_assessments()
+ assert ipp.results["Problem"][0]["Number of constraints"] == 61
+ assert ipp.results["Problem"][0]["Number of variables"] == 53
+ assert ipp.results["Problem"][0]["Number of nonzeros"] == 143
+
+ # *********************************************************************
+ # *********************************************************************
+
+ # validation
+
+ # only the IA arc should be installed
+ true_imp_arcs_selected = [True, False, False]
+ for node_pair, true_arc_decision in zip(list_imp_arcs, true_imp_arcs_selected):
+ assert (
+ true_arc_decision
+ in ipp.networks["mynet"]
+ .edges[(*node_pair, 0)][Network.KEY_ARC_TECH]
+ .options_selected
+ )
+ # only two arcs between A, B and C can be installed
+ arcs_selected = tuple(
+ 1
+ for node_pair in list_other_arcs
+ if True in ipp.networks["mynet"]
+ .edges[(*node_pair, 0)][Network.KEY_ARC_TECH]
+ .options_selected
+ )
+ assert sum(arcs_selected) == 2
+ # the network must be tree-shaped
+ assert ipp.networks["mynet"].has_tree_topology()
+ # capex
+ assert math.isclose(pyo.value(ipp.instance.var_capex), 10.0, abs_tol=1e-3)
+ # the objective function
+ assert math.isclose(pyo.value(ipp.instance.obj_f), -(10+(-11.93236715+10)), abs_tol=1e-3)
# *****************************************************************************
# *****************************************************************************
\ No newline at end of file
diff --git a/tests/test_gis_utils.py b/tests/test_gis_utils.py
index d2edd78..b561546 100644
--- a/tests/test_gis_utils.py
+++ b/tests/test_gis_utils.py
@@ -2650,15 +2650,35 @@ class TestGisUtils:
edge_key = (0, 1)
G.add_node(0, x=0, y=0)
G.add_node(1, x=1, y=1)
- G.add_edge(0, 1, length=2**0.5)
+ k1 = G.add_edge(0, 1, length=2**0.5)
first_geo = gis_utils.create_edge_geometry(G, edge_key)
edge_key = (2, 3)
G.add_node(2, x=0, y=1)
G.add_node(3, x=1, y=0)
- G.add_edge(2, 3, length=2**0.5)
+ k2 = G.add_edge(2, 3, length=2**0.5)
second_geo = gis_utils.create_edge_geometry(G, edge_key)
assert intersects(first_geo, second_geo)
+ # test finding overlapping edges
+ ol_edges = gis_utils.find_overlapping_edges(G)
+ assert len(ol_edges) == 1
+ assert ((0,1,k1),(2,3,k2)) in ol_edges
+ # no overlapping edges if there is only one node
+ ol_edges = gis_utils.find_overlapping_edges(G, excluded_edges=[(0,1,k1)])
+ assert len(ol_edges) == 0
+
+ # add fifth node
+ G.add_node(4, x=0.5, y=1)
+ k3 = G.add_edge(3, 4, length=(0.5**2+1**2)**0.5)
+ ol_edges = gis_utils.find_overlapping_edges(G)
+ assert len(ol_edges) == 2
+ assert ((0,1,k1),(2,3,k2)) in ol_edges
+ assert ((0,1,k1),(3,4,k3)) in ol_edges
+ # fewer overlapping edges
+ ol_edges = gis_utils.find_overlapping_edges(G, excluded_edges=[(2,3,k2)])
+ assert len(ol_edges) == 1
+ assert ((0,1,k1),(3,4,k3)) in ol_edges
+
# *************************************************************************
# *************************************************************************
--
GitLab