Fix several regressions in 'simplify()' - #166
Conversation
|
Hi @jnettels, sorry for messing this up. I probably rely to much on working tests, and was just content after my (green field) test case worked out well. I'll look into that immediately. I think, handling existing pipes (without investment) might be tackled the following way, you already partly implemented:
PS: A minimum working example for a case with existing pipes would be great. I tend to use it as an integration test if it's present. PPS: I extracted the problem you mentioned above from the dataset: roundtrip.zip PPPS: The problem occurs in the |
|
I think that solved one of the two remaining issues, nice! Now the "keeping existing pipes" To Do remains, will you tackle that, too? New test data: In the previous example input, I set two pipes as existing: This is the expected optimization result: The diagonal one in the middle has a large capacity, but just enough for all buildings connected directly. This forces usage of the outer path. There, the shorter detour is existing, but with a tiny capacity. Therefore the longer detour is needed for extra capacity To get this result, some of the changes in #164 are required. One last(?) thing we should do in here is drop the following in connect_points.py and remove # Keep only the shortest of all lines connecting the same two points
lines = go.drop_parallel_lines(lines)Your |
|
Thanks for the example. Reading through it made me reconsider my above suggestion: Treating existing pipes to have zero weight when searching for detours is only true if the existing capacity is sufficient. However, I do not see why existing pipes should be forced to stay in the optimisation. If it is cheaper to replace an existing pipe by a new one, why should you build a round-trip? So, I see the following cases:
Your suggestion, however is to always keep both branches if there is something existing on either side. This will delay the decision to a later stage and would be more flexible. I can live with that approach but would not force users to have existing pipes untouched in the optimisation. |
|
Hmm... interesting... to your points:
That assumption is the critical part. When you run the new example without the longer (non-existing) detour, you get the following result: This is true, because there simply is now way to provide capacity to the downstream consumers. This is the worst-case scenario from a "user experience" point of view, because
I would much rather have a solvable network. The user always has to check the results and then they will find the parallel paths. Then they have to make a decision of e.g. replacing the pipe in the input data.
For both cases you describe, the existing pipe may have a very large capacity (large DN). And the new detour might only need to be a small DN to achieve the total required capacity. There could always be cases where you would not declare the existing pipe "useless", if the investment for a complete replacement is much larger.
I must admit that I am not sure what you mean here. :-) |
Renaming back and forth is unnecessary.
|
Ok, then we agree. These "infeasible" errors were actually very annoying for me when I worked with a larger network where I had to manually check which consumers caused it. Now I realize that it should be possible to automatically create duplicates of the existing network segments and assign them a separate hp_type with much higher fixed costs. That could yield a result where the solver chooses them only if absolutely necessary. Then you could automatically determine where these existing and new pipes overlap and replace both with one pipe with their combined capacity, plus an additional attribute that flags a "replacement". ...well, that is a topic for some other time^^ |
This includes a refactor of _remove_useless_forks to allow reusing part of its logic.
|
I realised that the change to keep unique values can imply an issue: Now, there is a reason to have two edges connecting the same two nodes, which is currently not supported. |
|
I could hardly remember, but yes, that must have been the reason why I implemented However, here is what I found with our current test cases:
So, what is the consequence?
Do you think you can do that? |
In the process, I refactored some of the logic functions to be easier to understund. for example, `_have_unique_values` was True if there are different values. The new name inverts the logic to `_attribute_values_equal` and should be easier to understand. (Also, the separate NaN treatment in `_all_values_equal` was useless, as NaN != any value.)
|
Now, both subroutines of |
|
Hmmm.. we really need to step up those tests, yes. The "existing" detour now seems to be handled as required, nice! However, I found a new regression. The previous commit da0e64c was fine with Now we get
Sorry :-) PS A note to a previous comment:
Actually, in #164 one of my changes is storing the actual flow in the results. For existing pipes, a user can now compare capacity and flow. This allows them to directly see how much of the existing capacity is actually needed and identify pipe segments that are e.g. severely under-utilized. And yes, the simulation part is completely unsupported, to my knowledge. Plus pandapipes works well enough for what I would assume the simulation should have become. I have a bad feeling about deleting the code... but the tests should be disabled, I agree |
We do that in solph. But even if the tests are not ran automatically, you should really define (failing) tests. I ran all the example code you shared above with both datasets for street input data and they finished through just fine. Also, I do not get why
|
|
Ok, several things to discuss in parallel :-)
Agree! Once we are done and have finalized the "expected result", I would suggest adding the following to def test_process_geometry():
"""Test ``process_geometry()`` function with a simple example.
It includes pipes with ``existing=1``, due to which an otherwise
deleted detour must be kept.
"""
base_dir = os.path.join(os.path.dirname(__file__), "_files/process_geometry")
gdf_lines = gpd.read_file(os.path.join(base_dir, "in/lines_input_existing.geojson"))
gdf_prod = gpd.read_file(os.path.join(base_dir, "in/producers_polygon.geojson"))
gdf_cons = gpd.read_file(os.path.join(base_dir, "in/consumers_polygon.geojson"))
tn_input = cp.process_geometry(
lines=gdf_lines,
producers=gdf_prod,
consumers=gdf_cons,
method="boundary",
reset_index=True,
welding=True,
)
assert tn_input['pipes'].crs is not None
assert [c in tn_input['pipes'].columns for c in ['type', 'id_full']]
assert tn_input['pipes'].index.name == 'id'
# Update expected result
# tn_input['pipes'].to_file(os.path.join(base_dir, "out/pipes.geojson"))
# Load expected result
gdf_pipes_test = gpd.read_file(os.path.join(base_dir, "out/pipes.geojson")).set_index('id')
assert gdf_pipes_test.equals(tn_input['pipes'])
You are right, that is actually suspicious. I quickly did the following two changes:
For our simple test case here, our network solved correctly. We would have to ask Johannes why he implemented that restriction in the first place... but I do not see why the solver would care. There is a big HOWEVER, though, which leads to
The fact that the network can now be solved does not mean the result from Are you saying that you cannot reproduce the issue? import os
import io
import pandas as pd
import matplotlib.pyplot as plt
import geopandas as gpd
import dhnx
import logging
import dhnx
logging.basicConfig(
format='%(asctime)s %(name)-12s %(levelname)-8s %(message)s',
datefmt='%H:%M:%S')
logger = logging.getLogger(__name__) # Create a logger for this module
logger.setLevel(level='DEBUG') # Logger for this module
logging.getLogger('dhnx').setLevel(level='DEBUG')
def get_default_dhnx_invest_options():
"""Generate a dictionary with default investments options for DHNx."""
consumers_bus = """
label_2,active,excess,shortage,shortage costs,excess costs
heat,1,0,0,999999,9999
"""
consumers_demand = """
label_2,active,nominal_capacity
heat,1,1
"""
producers_bus = """
,label_2,active,excess,shortage,shortage costs,excess costs
1,heat,1,0,0,9999,9999
"""
producers_source = """
label_2,active
heat,1
"""
network = """
label_3,nonconvex,l_factor,l_factor_fix,cap_min,cap_max,capex_pipes,fix_costs
pipe-generic,1,1.5812551e-07,0.0221164,1,519307.085403,0.02908927,948.69105059
"""
invest_options = dict(
consumers=dict(
bus=pd.read_csv(io.StringIO(consumers_bus)),
demand=pd.read_csv(io.StringIO(consumers_demand)),
),
producers=dict(
bus=pd.read_csv(io.StringIO(producers_bus)),
source=pd.read_csv(io.StringIO(producers_source)),
),
network=dict(
pipes=pd.read_csv(io.StringIO(network))
),
)
return invest_options
# Load input
gdf_lines_streets = gpd.read_file("./input/streets_input.geojson")
# gdf_lines_streets = gpd.read_file("./input/streets_input_existing.geojson")
# gdf_lines_streets = gpd.read_file("./input/streets_input_detour_merged.geojson")
gdf_poly_gen = gpd.read_file("./input/producers_polygon.geojson")
gdf_poly_houses = gpd.read_file("./input/consumers_polygon.geojson")
# Plot input
_, ax = plt.subplots(figsize=(20, 10), dpi=300)
gdf_poly_houses.to_crs("EPSG:4647").plot(ax=ax, color='green')
gdf_poly_gen.to_crs("EPSG:4647").plot(ax=ax, color='red')
gdf_lines_streets.to_crs("EPSG:4647").plot(ax=ax, color='blue')
plt.title('Input geometry')
plt.show()
gdf_lines_streets = gdf_lines_streets.replace({'hp_type': {'pipe-typ-A': 'pipe-generic'}})
gdf_lines_streets = gdf_lines_streets.fillna({'existing': 0})
gdf_lines_streets = gdf_lines_streets.replace({None: float('nan')})
if 'existing' in gdf_lines_streets.columns:
_, ax = plt.subplots(figsize=(20, 10), dpi=300)
gdf_lines_streets.plot(ax=ax, column='existing', legend_kwds={'label': 'existing'}, legend=True)
plt.title('Status "existing" of input pipes')
plt.show()
# Process the geometry
reset_index = True
tn_input = dhnx.gistools.connect_points.process_geometry(
lines=gdf_lines_streets,
producers=gdf_poly_gen.copy(),
consumers=gdf_poly_houses.copy(),
method="boundary",
reset_index=reset_index,
# welding=False,
)
# Plot output after processing the geometry
_, ax = plt.subplots(figsize=(20, 10), dpi=300)
tn_input['consumers'].plot(ax=ax, color='green')
tn_input['producers'].plot(ax=ax, color='red')
tn_input['pipes'].plot(ax=ax, color='blue')
tn_input['forks'].plot(ax=ax, color='grey')
plt.title('Geometry after pre-processing')
plt.show()
# breakpoint()
print("Tests:")
print(tn_input['pipes'].crs is not None)
print([c in tn_input['pipes'].columns for c in ['type', 'id_full']])
if reset_index:
print(tn_input['pipes'].index.name == 'id')
# Optionally export the geodataframes and load it into qgis
# for checking the results of the geometry processing
save_path = './output'
if save_path is not None:
if not os.path.exists(os.path.dirname(save_path)):
os.makedirs(os.path.dirname(save_path))
for filename, gdf in tn_input.items():
try:
gdf.to_file(os.path.join(save_path, filename)+'.geojson')
gdf.to_file(os.path.join(save_path, filename)+'.gpkg')
except Exception as e:
print(gdf)
breakpoint()
logger.exception(e)
# initialize a ThermalNetwork
network = dhnx.network.ThermalNetwork()
# add the pipes, forks, consumer, and producers to the ThermalNetwork
for k, v in tn_input.items():
network.components[k] = v
# check if ThermalNetwork is consistent
network.is_consistent()
# load the specification of the oemof-solph components
invest_opt = get_default_dhnx_invest_options()
settings = dict(
# solver='cbc',
solver='gurobi',
# return_existing=True, # Requires PR # 164
)
# Perform the investment optimisation
network.optimize_investment(invest_options=invest_opt, **settings)
results_edges = network.results.optimization['components']['pipes']
gdf_pipes = network.components['pipes'].copy()
cols_drop = [c for c in results_edges.columns if c in gdf_pipes]
gdf_pipes = gdf_pipes.drop(columns=cols_drop) # Drop duplicate columns
gdf_pipes = gdf_pipes.join(results_edges, rsuffix='_results')
gdf_pipes = gdf_pipes[gdf_pipes['capacity'] > 0] # Keep only DN>0 in output
# Plot output after processing the geometry
_, ax = plt.subplots(figsize=(20, 10), dpi=300)
gdf_poly_houses.to_crs("EPSG:4647").plot(ax=ax, color='green')
gdf_poly_gen.to_crs("EPSG:4647").plot(ax=ax, color='red')
gdf_pipes.to_crs("EPSG:4647").plot(ax=ax, column='capacity')
plt.title('Network result')
plt.show()
gdf_pipes.to_file('./output/pipes_result.gpkg')
_, ax = plt.subplots(figsize=(20, 10), dpi=300)
gdf_pipes.to_crs("EPSG:4647").plot(ax=ax, column='type', legend=True, cmap='viridis')
gdf_poly_houses.to_crs("EPSG:4647").plot(ax=ax, color='green')
gdf_poly_gen.to_crs("EPSG:4647").plot(ax=ax, color='red')
plt.title('Line type (distribution/generator/house)')
plt.show() |
It allows construction of a "reverse" house connection lines, and was provided in #166 (comment)
There were two problems when inverting "_all_values_equal": First, the logic in _remove_useless_forks was notinverted. This is now more transparent without continue statement. Second, series of all NaN were not merged because NaN != NaN.
Coverage reportClick to see where and how coverage changed
This report was generated by python-coverage-comment-action |
||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
Turns out, what I considered to be refactoring actually broke @jnettels: This should be done, now. Just one thing might still need discussion:
This is true, but simplification is not mandatory. However, the old functionality was also broken (in the sense of the original |
This is not needed anymore because simplify() already includes this step: It uses a networkx Graph, which ony allows one edge between two nodes.
|
Awesome! Thanks for getting all the tests running!
Agreed I hope I do not mess up this branch by also pushing to it, but I will risk it:
In 236505f I added my updated test case as a new test for Maybe this is the solution: Once we get an expected result, we need to comment out How do you feel about renaming the |
Indeed this seems to be the only change required to fix the error and get an expected output. |
|
Hmm... the actual error in But my tests fail when comparing with the expected |
|
Are you aware that |
|
Good point. But The failing columns were The only remaining question from me would be: |
|
Only doctest link check fails because of traffic limitations. I will merge. |







Hi @p-snft, now that I have had an opportunity to review #157, these are the possible improvements and regressions that I found. I consider this a blocker that needs to be solved before I can continue with #164. I considered putting all the changes in there but I hope it is cleaner this way.
TLDR: I would like to assign the remaining work back to you.
Stuff already fixed
simplify()pathin the gdf returned bysimplify(). It adds unnescessary size to the GeoDataFramesimplify()compared to the oldweld_segments():linesis lost in the processlinesis losttypeandid_fullof GeoDataFramelineswere just dropped (which is related to the next point)dhnx to fail, because it expects connection lines to always go towards producers/consumers.
(see e.g. dhs_nodes.py line 170
"Pipes must not go from 'consumers'!")and building connection lines.
"Luckily", it is the same problem I need to solve to enable support for existing pipes.
For given columns (here
type), pipes must not be merged if their values ('DL', 'GL', 'HL') differ.Most of the above were quite simple fixes. Only the function to prevent certain pipes from merging was a little bit more involved. I hope my solution with
retain_unique_valuesis acceptable.Remaining issues
weldingtosimplifygeometry_operations.drop_parallel_lines(), as it is superseded bysimplify()existing=1pipe (i.e. a pipe with a fixed capacity). We cannot know up front if the path containing the existing pipe has the required capacity to allow a feasible network. When working with existing pipes you sometimes have to deliberately create detours that the solver will use once the capacities of the existing pipes are exhausted.simplifyjust does not seem to do what it is supposed to, i.e. of two alternatives a longer path is preserved. See the example belowConsider this input:

Commit 1e2b011 (last version with old 'welding', before merge of new simplify()) yields

This is fine, only the roundabout could potentially be simplified
The result of the current dev branch 465d067 and after including my fixes is the same:

I messed with the input geometry a bit, removing and redrawing that detour by hand, but the results were the same.
At this point I would like to assign this back to you, if possible. I have attached the python script and input data for reproducing the issue. I hope you have a better chance to fix it.
I hope the download works:
test_dhnx_simplify_regressions.zip