s_ranjbar/energy_system_modelling_workflow
2
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

feature: add district heating network sizing workflow

Browse Source 
Reviewed-on: https://nextgenerations-cities.encs.concordia.ca/gitea/s_ranjbar/energy_system_modelling_workflow/pulls/18
This commit is contained in:
Majid Rezaei 2024-08-15 15:50:16 -04:00
commit 9af87fe482
10 changed files with 1586035 additions and 73 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

111
district_heating_network.py Normal file
View File

@ -0,0 +1,111 @@
from scripts.district_heating_network.directory_manager import DirectoryManager
import subprocess
from scripts.ep_run_enrich import energy_plus_workflow
from hub.imports.geometry_factory import GeometryFactory
from hub.helpers.dictionaries import Dictionaries
from hub.imports.construction_factory import ConstructionFactory
from hub.imports.usage_factory import UsageFactory
from hub.imports.weather_factory import WeatherFactory
from hub.imports.results_factory import ResultFactory
from scripts.energy_system_retrofit_report import EnergySystemRetrofitReport
from scripts.geojson_creator import process_geojson
from scripts import random_assignation
from hub.imports.energy_systems_factory import EnergySystemsFactory
from scripts.energy_system_sizing import SystemSizing
from scripts.solar_angles import CitySolarAngles
from scripts.pv_sizing_and_simulation import PVSizingSimulation
from scripts.energy_system_retrofit_results import consumption_data, cost_data
from scripts.energy_system_sizing_and_simulation_factory import EnergySystemsSimulationFactory
from scripts.costs.cost import Cost
from scripts.costs.constants import SKIN_RETROFIT_AND_SYSTEM_RETROFIT_AND_PV, SYSTEM_RETROFIT_AND_PV, CURRENT_STATUS
import hub.helpers.constants as cte
from hub.exports.exports_factory import ExportsFactory
from scripts.pv_feasibility import pv_feasibility
import matplotlib.pyplot as plt
import numpy as np
from scripts.district_heating_network.district_heating_network_creator import DistrictHeatingNetworkCreator
from scripts.district_heating_network.district_heating_factory import DistrictHeatingFactory
import json
#%% --------------------------------------------------------------------------------------------------------------------
# Manage File Path
base_path = "./"
dir_manager = DirectoryManager(base_path)
# Input files directory
input_files_path = dir_manager.create_directory('input_files')
geojson_file_path = input_files_path / 'output_buildings.geojson'
pipe_data_file = input_files_path / 'pipe_data.json'
# Output files directory
output_path = dir_manager.create_directory('out_files')
# Subdirectories for output files
energy_plus_output_path = dir_manager.create_directory('out_files/energy_plus_outputs')
simulation_results_path = dir_manager.create_directory('out_files/simulation_results')
sra_output_path = dir_manager.create_directory('out_files/sra_outputs')
cost_analysis_output_path = dir_manager.create_directory('out_files/cost_analysis')
#%% --------------------------------------------------------------------------------------------------------------------
# Area Under Study
location = [45.4934614681437, -73.57982834742518]
#%% --------------------------------------------------------------------------------------------------------------------
# Create geojson of buildings
process_geojson(x=location[1], y=location[0], diff=0.001)
#%% --------------------------------------------------------------------------------------------------------------------
# Create ciry and run energyplus workflow
city = GeometryFactory(file_type='geojson',
path=geojson_file_path,
height_field='height',
year_of_construction_field='year_of_construction',
function_field='function',
function_to_hub=Dictionaries().montreal_function_to_hub_function).city
ConstructionFactory('nrcan', city).enrich()
UsageFactory('nrcan', city).enrich()
WeatherFactory('epw', city).enrich()
# SRA
ExportsFactory('sra', city, output_path).export()
sra_path = (output_path / f'{city.name}_sra.xml').resolve()
subprocess.run(['sra', str(sra_path)])
ResultFactory('sra', city, output_path).enrich()
# EP Workflow
energy_plus_workflow(city, energy_plus_output_path)
#%% --------------------------------------------------------------------------------------------------------------------
# District Heating Network Creator
central_plant_locations = [(-73.57812571080625, 45.49499447346277)] # Add at least one location
roads_file = "./input_files/roads.json"
dhn_creator = DistrictHeatingNetworkCreator(geojson_file_path, roads_file, central_plant_locations)
network_graph = dhn_creator.run()
#%% --------------------------------------------------------------------------------------------------------------------
# Pipe and pump sizing
with open(pipe_data_file, 'r') as f:
pipe_data = json.load(f)
factory = DistrictHeatingFactory(
city=city,
graph=network_graph,
supply_temperature=80 + 273, # in Kelvin
return_temperature=60 + 273, # in Kelvin
simultaneity_factor=0.9
)
factory.enrich()
factory.sizing()
factory.calculate_diameters_and_costs(pipe_data)
pipe_groups, total_cost = factory.analyze_costs()
# Save the pipe groups with total costs to a CSV file
factory.save_pipe_groups_to_csv('pipe_groups.csv')
#%% --------------------------------------------------------------------------------------------------------------------

191
input_files/pipe_data.json Normal file
View File

@ -0,0 +1,191 @@
[
{
"DN": 16,
"inner_diameter": 16.1,
"outer_diameter": 21.3,
"thickness": 2.6,
"cost_per_meter": 320
},
{
"DN": 20,
"inner_diameter": 21.7,
"outer_diameter": 26.9,
"thickness": 2.6,
"cost_per_meter": 320
},
{
"DN": 25,
"inner_diameter": 27.3,
"outer_diameter": 33.7,
"thickness": 3.2,
"cost_per_meter": 320
},
{
"DN": 32,
"inner_diameter": 37.2,
"outer_diameter": 42.4,
"thickness": 2.6,
"cost_per_meter": 350
},
{
"DN": 40,
"inner_diameter": 43.1,
"outer_diameter": 48.3,
"thickness": 2.6,
"cost_per_meter": 375
},
{
"DN": 50,
"inner_diameter": 54.5,
"outer_diameter": 60.3,
"thickness": 2.9,
"cost_per_meter": 400
},
{
"DN": 65,
"inner_diameter": 70.3,
"outer_diameter": 76.1,
"thickness": 2.9,
"cost_per_meter": 450
},
{
"DN": 80,
"inner_diameter": 82.5,
"outer_diameter": 88.9,
"thickness": 3.2,
"cost_per_meter": 480
},
{
"DN": 90,
"inner_diameter": 100.8,
"outer_diameter": 108,
"thickness": 3.6,
"cost_per_meter": 480
},
{
"DN": 100,
"inner_diameter": 107.1,
"outer_diameter": 114.3,
"thickness": 3.6,
"cost_per_meter": 550
},
{
"DN": 110,
"inner_diameter": 125.8,
"outer_diameter": 133,
"thickness": 3.6,
"cost_per_meter": 550
},
{
"DN": 125,
"inner_diameter": 132.5,
"outer_diameter": 139.7,
"thickness": 3.6,
"cost_per_meter": 630
},
{
"DN": 140,
"inner_diameter": 151,
"outer_diameter": 159,
"thickness": 4,
"cost_per_meter": 700
},
{
"DN": 150,
"inner_diameter": 160.3,
"outer_diameter": 168.3,
"thickness": 4,
"cost_per_meter": 700
},
{
"DN": 180,
"inner_diameter": 184.7,
"outer_diameter": 193.7,
"thickness": 4.5,
"cost_per_meter": 700
},
{
"DN": 200,
"inner_diameter": 210.1,
"outer_diameter": 219.1,
"thickness": 4.5,
"cost_per_meter": 860
},
{
"DN": 250,
"inner_diameter": 263,
"outer_diameter": 273,
"thickness": 5,
"cost_per_meter": 860
},
{
"DN": 300,
"inner_diameter": 312.7,
"outer_diameter": 323.9,
"thickness": 5.6,
"cost_per_meter": 860
},
{
"DN": 350,
"inner_diameter": 344.4,
"outer_diameter": 355.6,
"thickness": 5.6,
"cost_per_meter": 860
},
{
"DN": 400,
"inner_diameter": 393.8,
"outer_diameter": 406.4,
"thickness": 6.3,
"cost_per_meter": 860
},
{
"DN": 450,
"inner_diameter": 444.4,
"outer_diameter": 457,
"thickness": 6.3,
"cost_per_meter": 860
},
{
"DN": 500,
"inner_diameter": 495.4,
"outer_diameter": 508,
"thickness": 6.3,
"cost_per_meter": 860
},
{
"DN": 600,
"inner_diameter": 595.8,
"outer_diameter": 610,
"thickness": 7.1,
"cost_per_meter": 860
},
{
"DN": 700,
"inner_diameter": 696.8,
"outer_diameter": 711,
"thickness": 7.1,
"cost_per_meter": 860
},
{
"DN": 800,
"inner_diameter": 797,
"outer_diameter": 813,
"thickness": 8,
"cost_per_meter": 860
},
{
"DN": 900,
"inner_diameter": 894,
"outer_diameter": 914,
"thickness": 10,
"cost_per_meter": 860
},
{
"DN": 1000,
"inner_diameter": 996,
"outer_diameter": 1016,
"thickness": 10,
"cost_per_meter": 860
}
]

1585013
input_files/roads.json Normal file
View File

File diff suppressed because it is too large Load Diff

43
main.py
View File

@ -1,4 +1,5 @@
from pathlib import Path
from scripts.district_heating_network.directory_manager import DirectoryManager
import subprocess
from scripts.ep_run_enrich import energy_plus_workflow
from hub.imports.geometry_factory import GeometryFactory
@ -22,23 +23,31 @@ import hub.helpers.constants as cte
from hub.exports.exports_factory import ExportsFactory
from scripts.pv_feasibility import pv_feasibility
import matplotlib.pyplot as plt
import numpy as np
# Specify the GeoJSON file path
data = {}
input_files_path = (Path(__file__).parent / 'input_files')
input_files_path.mkdir(parents=True, exist_ok=True)
geojson_file = process_geojson(x=-73.58001358793511, y=45.496445294438715, diff=0.0001)
from scripts.district_heating_network.district_heating_network_creator import DistrictHeatingNetworkCreator
from scripts.district_heating_network.road_processor import road_processor
from scripts.district_heating_network.district_heating_factory import DistrictHeatingFactory
base_path = Path(__file__).parent
dir_manager = DirectoryManager(base_path)
# Input files directory
input_files_path = dir_manager.create_directory('input_files')
geojson_file_path = input_files_path / 'output_buildings.geojson'
output_path = (Path(__file__).parent / 'out_files').resolve()
output_path.mkdir(parents=True, exist_ok=True)
energy_plus_output_path = output_path / 'energy_plus_outputs'
energy_plus_output_path.mkdir(parents=True, exist_ok=True)
simulation_results_path = (Path(__file__).parent / 'out_files' / 'simulation_results').resolve()
simulation_results_path.mkdir(parents=True, exist_ok=True)
sra_output_path = output_path / 'sra_outputs'
sra_output_path.mkdir(parents=True, exist_ok=True)
cost_analysis_output_path = output_path / 'cost_analysis'
cost_analysis_output_path.mkdir(parents=True, exist_ok=True)
# Output files directory
output_path = dir_manager.create_directory('out_files')
# Subdirectories for output files
energy_plus_output_path = dir_manager.create_directory('out_files/energy_plus_outputs')
simulation_results_path = dir_manager.create_directory('out_files/simulation_results')
sra_output_path = dir_manager.create_directory('out_files/sra_outputs')
cost_analysis_output_path = dir_manager.create_directory('out_files/cost_analysis')
# Select city area
location = [45.53067276979674, -73.70234652694087]
process_geojson(x=location[1], y=location[0], diff=0.001)
# Create city object
city = GeometryFactory(file_type='geojson',
path=geojson_file_path,
height_field='height',
@ -83,4 +92,4 @@ UsageFactory('nrcan', city).enrich()
# # Save the plot
# plt.savefig('plot_nrcan.png', bbox_inches='tight')
# plt.close()
print('test')
print('test')

4
pv_assessment.py
View File

@ -14,7 +14,7 @@ import hub.helpers.constants as cte
from hub.exports.exports_factory import ExportsFactory
from scripts.pv_sizing_and_simulation import PVSizingSimulation
# Specify the GeoJSON file path
geojson_file = process_geojson(x=-73.5681295982132, y=45.49218262677643, diff=0.0005)
geojson_file = process_geojson(x=-73.5681295982132, y=45.49218262677643, diff=0.0001)
file_path = (Path(__file__).parent / 'input_files' / 'output_buildings.geojson')
# Specify the output path for the PDF file
output_path = (Path(__file__).parent / 'out_files').resolve()
@ -34,7 +34,7 @@ ExportsFactory('sra', city, output_path).export()
sra_path = (output_path / f'{city.name}_sra.xml').resolve()
subprocess.run(['sra', str(sra_path)])
ResultFactory('sra', city, output_path).enrich()
energy_plus_workflow(city)
energy_plus_workflow(city, output_path=output_path)
solar_angles = CitySolarAngles(city.name,
city.latitude,
city.longitude,

16
scripts/district_heating_network/directory_manager.py Normal file
View File

@ -0,0 +1,16 @@
from pathlib import Path
class DirectoryManager:
def __init__(self, base_path):
self.base_path = Path(base_path)
self.directories = {}
def create_directory(self, relative_path):
full_path = self.base_path / relative_path
full_path.mkdir(parents=True, exist_ok=True)
self.directories[relative_path] = full_path
return full_path
def get_directory(self, relative_path):
return self.directories.get(relative_path, None)

248
scripts/district_heating_network/district_heating_factory.py Normal file
View File

@ -0,0 +1,248 @@
import CoolProp.CoolProp as CP
import math
import logging
import numpy as np
import csv
class DistrictHeatingFactory:
"""
DistrictHeatingFactory class
This class is responsible for managing the district heating network, including
enriching the network graph with building data, calculating flow rates,
sizing pipes, and analyzing costs.
"""
def __init__(self, city, graph, supply_temperature, return_temperature, simultaneity_factor):
"""
Initialize the DistrictHeatingFactory object.
:param city: The city object containing buildings and their heating demands.
:param graph: The network graph representing the district heating network.
:param supply_temperature: The supply temperature of the heating fluid in the network (°C).
:param return_temperature: The return temperature of the heating fluid in the network (°C).
:param simultaneity_factor: The simultaneity factor used to adjust flow rates for non-building pipes.
"""
self._city = city
self._network_graph = graph
self._supply_temperature = supply_temperature
self._return_temperature = return_temperature
self.simultaneity_factor = simultaneity_factor
self.fluid = "Water" # The fluid used in the heating network
def enrich(self):
"""
Enrich the network graph nodes with the whole building object from the city buildings.
This method associates each building node in the network graph with its corresponding
building object from the city, allowing access to heating demand data during calculations.
"""
for node_id, node_attrs in self._network_graph.nodes(data=True):
if node_attrs.get('type') == 'building':
building_name = node_attrs.get('name')
building_found = False
for building in self._city.buildings:
if building.name == building_name:
self._network_graph.nodes[node_id]['building_obj'] = building
building_found = True
break
if not building_found:
logging.error(msg=f"Building with name '{building_name}' not found in city.")
def calculate_flow_rates(self, A, Gext):
"""
Solve the linear system to find the flow rates in each branch.
:param A: The incidence matrix representing the network connections.
:param Gext: The external flow rates for each node in the network.
:return: The calculated flow rates for each edge, or None if an error occurs.
"""
try:
G = np.linalg.lstsq(A, Gext, rcond=None)[0]
return G
except np.linalg.LinAlgError as e:
logging.error(f"Error solving the linear system: {e}")
return None
def switch_nodes(self, A, edge_index, node_index, edge):
"""
Switch the in and out nodes for the given edge in the incidence matrix A.
:param A: The incidence matrix representing the network connections.
:param edge_index: The index of edges in the incidence matrix.
:param node_index: The index of nodes in the incidence matrix.
:param edge: The edge (u, v) to switch.
"""
u, v = edge
i = node_index[u]
j = node_index[v]
k = edge_index[edge]
A[i, k], A[j, k] = -A[i, k], -A[j, k]
def sizing(self):
"""
Calculate the hourly mass flow rates, assign them to the edges, and determine the pipe diameters.
This method generates the flow rates for each hour, adjusting the incidence matrix as needed to
ensure all flow rates are positive. It also applies the simultaneity factor to non-building pipes.
"""
num_nodes = self._network_graph.number_of_nodes()
num_edges = self._network_graph.number_of_edges()
A = np.zeros((num_nodes, num_edges)) # Initialize incidence matrix
node_index = {node: i for i, node in enumerate(self._network_graph.nodes())}
edge_index = {edge: i for i, edge in enumerate(self._network_graph.edges())}
# Initialize mass flow rate attribute for each edge
for u, v, data in self._network_graph.edges(data=True):
self._network_graph.edges[u, v]['mass_flow_rate'] = {"hour": [], "peak": None}
# Get the length of the hourly demand for the first building (assuming all buildings have the same length)
building = next(iter(self._city.buildings))
num_hours = len(building.heating_demand['hour'])
# Loop through each hour to generate Gext and solve AG = Gext
for hour in range(8760):
Gext = np.zeros(num_nodes)
# Calculate the hourly mass flow rates for each edge and fill Gext
for edge in self._network_graph.edges(data=True):
u, v, data = edge
for node in [u, v]:
if self._network_graph.nodes[node].get('type') == 'building':
building = self._network_graph.nodes[node].get('building_obj')
if building and "hour" in building.heating_demand:
hourly_demand = building.heating_demand["hour"][hour] # Get demand for current hour
specific_heat_capacity = CP.PropsSI('C', 'T', (self._supply_temperature + self._return_temperature) / 2,
'P', 101325, self.fluid)
mass_flow_rate = hourly_demand / 3600 / (
specific_heat_capacity * (self._supply_temperature - self._return_temperature))
Gext[node_index[node]] += mass_flow_rate
# Update incidence matrix A
i = node_index[u]
j = node_index[v]
k = edge_index[(u, v)]
A[i, k] = 1
A[j, k] = -1
# Solve for G (flow rates)
G = self.calculate_flow_rates(A, Gext)
if G is None:
return
# Check for negative flow rates and adjust A accordingly
iterations = 0
max_iterations = num_edges * 2
while any(flow_rate < 0 for flow_rate in G) and iterations < max_iterations:
for idx, flow_rate in enumerate(G):
if flow_rate < 0:
G[idx] = -G[idx] # Invert the sign directly
iterations += 1
# Store the final flow rates in the edges for this hour
for idx, (edge, flow_rate) in enumerate(zip(self._network_graph.edges(), G)):
u, v = edge
if not (self._network_graph.nodes[u].get('type') == 'building' or self._network_graph.nodes[v].get(
'type') == 'building'):
flow_rate *= self.simultaneity_factor # Apply simultaneity factor for non-building pipes
data = self._network_graph.edges[u, v]
data['mass_flow_rate']["hour"].append(flow_rate) # Append the calculated flow rate
# Calculate the peak flow rate for each edge
for u, v, data in self._network_graph.edges(data=True):
data['mass_flow_rate']['peak'] = max(data['mass_flow_rate']['hour'])
def calculate_diameters_and_costs(self, pipe_data):
"""
Calculate the diameter and costs of the pipes based on the maximum flow rate in each edge.
:param pipe_data: A list of dictionaries containing pipe specifications, including inner diameters
and costs per meter for different nominal diameters (DN).
"""
for u, v, data in self._network_graph.edges(data=True):
flow_rate = data.get('mass_flow_rate', {}).get('peak')
if flow_rate is not None:
try:
# Calculate the density of the fluid
density = CP.PropsSI('D', 'T', (self._supply_temperature + self._return_temperature) / 2, 'P', 101325,
self.fluid)
velocity = 0.9 # Desired fluid velocity in m/s
# Calculate the diameter of the pipe required for the given flow rate
diameter = math.sqrt((4 * abs(flow_rate)) / (density * velocity * math.pi)) * 1000 # Convert to mm
self._network_graph.edges[u, v]['diameter'] = diameter
# Match to the closest nominal diameter from the pipe data
closest_pipe = self.match_nominal_diameter(diameter, pipe_data)
self._network_graph.edges[u, v]['nominal_diameter'] = closest_pipe['DN']
self._network_graph.edges[u, v]['cost_per_meter'] = closest_pipe['cost_per_meter']
except Exception as e:
logging.error(f"Error calculating diameter or matching nominal diameter for edge ({u}, {v}): {e}")
def match_nominal_diameter(self, diameter, pipe_data):
"""
Match the calculated diameter to the closest nominal diameter.
:param diameter: The calculated diameter of the pipe (in mm).
:param pipe_data: A list of dictionaries containing pipe specifications, including inner diameters
and costs per meter for different nominal diameters (DN).
:return: The dictionary representing the pipe with the closest nominal diameter.
"""
closest_pipe = min(pipe_data, key=lambda x: abs(x['inner_diameter'] - diameter))
return closest_pipe
def analyze_costs(self):
"""
Analyze the costs based on the nominal diameters of the pipes.
This method calculates the total cost of piping for each nominal diameter group
and returns a summary of the grouped pipes and the total cost.
:return: A tuple containing the grouped pipe data and the total cost of piping.
"""
pipe_groups = {}
total_cost = 0 # Initialize total cost
for u, v, data in self._network_graph.edges(data=True):
dn = data.get('nominal_diameter')
if dn is not None:
pipe_length = self._network_graph.edges[u, v].get('length', 1) * 2 # Multiply by 2 for supply and return
cost_per_meter = data.get('cost_per_meter', 0)
if dn not in pipe_groups:
pipe_groups[dn] = {
'DN': dn,
'total_length': 0,
'cost_per_meter': cost_per_meter
}
pipe_groups[dn]['total_length'] += pipe_length
group_cost = pipe_length * cost_per_meter
total_cost += group_cost # Add to total cost
# Calculate total cost for each group
for group in pipe_groups.values():
group['total_cost'] = group['total_length'] * group['cost_per_meter']
return pipe_groups, total_cost # Return both the grouped data and total cost
def save_pipe_groups_to_csv(self, filename):
"""
Save the pipe groups and their total lengths to a CSV file.
:param filename: The name of the CSV file to save the data to.
"""
pipe_groups, _ = self.analyze_costs()
with open(filename, mode='w', newline='') as file:
writer = csv.writer(file)
# Write the header
writer.writerow(["Nominal Diameter (DN)", "Total Length (m)", "Cost per Meter", "Total Cost"])
# Write the data for each pipe group
for group in pipe_groups.values():
writer.writerow([
group['DN'],
group['total_length'],
group['cost_per_meter'],
group['total_cost']
])

372
scripts/district_heating_network/district_heating_network_creator.py Normal file
View File

@ -0,0 +1,372 @@
import json
import math
import logging
import matplotlib.pyplot as plt
import networkx as nx
from shapely.geometry import Polygon, Point, LineString
from typing import List, Tuple
from rtree import index
# Configure logging
logging.basicConfig(level=logging.INFO, format="%(asctime)s - %(levelname)s - %(message)s")
logging.getLogger("numexpr").setLevel(logging.ERROR)
def haversine(lon1, lat1, lon2, lat2):
"""
Calculate the great-circle distance between two points
on the Earth specified by their longitude and latitude.
"""
R = 6371000 # Radius of the Earth in meters
phi1 = math.radians(lat1)
phi2 = math.radians(lat2)
delta_phi = math.radians(lat2 - lat1)
delta_lambda = math.radians(lon2 - lon1)
a = math.sin(delta_phi / 2.0) ** 2 + \
math.cos(phi1) * math.cos(phi2) * \
math.sin(delta_lambda / 2.0) ** 2
c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))
return R * c # Output distance in meters
class DistrictHeatingNetworkCreator:
def __init__(self, buildings_file: str, roads_file: str, central_plant_locations: List[Tuple[float, float]]):
"""
Initialize the class with paths to the buildings and roads data files, and central plant locations.
:param buildings_file: Path to the GeoJSON file containing building data.
:param roads_file: Path to the GeoJSON file containing roads data.
:param central_plant_locations: List of tuples containing the coordinates of central plant locations.
"""
if len(central_plant_locations) < 1:
raise ValueError("The list of central plant locations must have at least one member.")
self.buildings_file = buildings_file
self.roads_file = roads_file
self.central_plant_locations = central_plant_locations
def run(self) -> nx.Graph:
"""
Main method to execute the district heating network creation process.
:return: NetworkX graph with nodes and edges representing the network.
"""
try:
self._load_and_process_data()
self._find_nearest_roads()
self._find_nearest_points()
self._break_down_roads()
self._create_graph()
self._create_mst()
self._iteratively_remove_edges()
self._add_centroids_to_mst()
self._convert_edge_weights_to_meters()
self._create_final_network_graph()
return self.network_graph
except Exception as e:
logging.error(f"Error during network creation: {e}")
raise
def _load_and_process_data(self):
"""
Load and process the building and road data.
"""
try:
# Load building data
with open(self.buildings_file, 'r') as file:
city = json.load(file)
self.centroids = []
self.building_names = []
self.building_positions = []
buildings = city['features']
for building in buildings:
coordinates = building['geometry']['coordinates'][0]
building_polygon = Polygon(coordinates)
centroid = building_polygon.centroid
self.centroids.append(centroid)
self.building_names.append(str(building['id']))
self.building_positions.append((centroid.x, centroid.y))
# Add central plant locations as centroids
for i, loc in enumerate(self.central_plant_locations, start=1):
centroid = Point(loc)
self.centroids.append(centroid)
self.building_names.append(f'central_plant_{i}')
self.building_positions.append((centroid.x, centroid.y))
# Load road data
with open(self.roads_file, 'r') as file:
roads = json.load(file)
line_features = [feature for feature in roads['features'] if feature['geometry']['type'] == 'LineString']
self.lines = [LineString(feature['geometry']['coordinates']) for feature in line_features]
self.cleaned_lines = [LineString([line.coords[0], line.coords[-1]]) for line in self.lines]
except Exception as e:
logging.error(f"Error loading and processing data: {e}")
raise
def _find_nearest_roads(self):
"""
Find the nearest road for each building centroid.
"""
try:
self.closest_roads = []
unique_roads_set = set()
# Create spatial index for roads
idx = index.Index()
for pos, line in enumerate(self.cleaned_lines):
idx.insert(pos, line.bounds)
for centroid in self.centroids:
min_distance = float('inf')
closest_road = None
for pos in idx.nearest(centroid.bounds, 10):
road = self.cleaned_lines[pos]
distance = road.distance(centroid)
if distance < min_distance:
min_distance = distance
closest_road = road
if closest_road and closest_road.wkt not in unique_roads_set:
unique_roads_set.add(closest_road.wkt)
self.closest_roads.append(closest_road)
except Exception as e:
logging.error(f"Error finding nearest roads: {e}")
raise
def _find_nearest_points(self):
"""
Find the nearest point on each closest road for each centroid.
"""
def find_nearest_point_on_line(point: Point, line: LineString) -> Point:
return line.interpolate(line.project(point))
try:
self.nearest_points = []
for centroid in self.centroids:
min_distance = float('inf')
closest_road = None
for road in self.closest_roads:
distance = centroid.distance(road)
if distance < min_distance:
min_distance = distance
closest_road = road
if closest_road:
nearest_point = find_nearest_point_on_line(centroid, closest_road)
self.nearest_points.append(nearest_point)
except Exception as e:
logging.error(f"Error finding nearest points: {e}")
raise
def _break_down_roads(self):
"""
Break down roads into segments connecting nearest points.
"""
def break_down_roads(closest_roads: List[LineString], nearest_points_list: List[Point]) -> List[LineString]:
new_segments = []
for road in closest_roads:
coords = list(road.coords)
points_on_road = [point for point in nearest_points_list if road.distance(point) < 0.000000001]
sorted_points = sorted(points_on_road, key=lambda point: road.project(point))
sorted_points.insert(0, Point(coords[0]))
sorted_points.append(Point(coords[-1]))
for i in range(len(sorted_points) - 1):
segment = LineString([sorted_points[i], sorted_points[i + 1]])
new_segments.append(segment)
return new_segments
try:
self.new_segments = break_down_roads(self.closest_roads, self.nearest_points)
self.cleaned_lines = [line for line in self.cleaned_lines if line not in self.closest_roads]
self.cleaned_lines.extend(self.new_segments)
except Exception as e:
logging.error(f"Error breaking down roads: {e}")
raise
def _create_graph(self):
"""
Create a NetworkX graph from the cleaned lines.
"""
try:
self.G = nx.Graph()
for line in self.cleaned_lines:
coords = list(line.coords)
for i in range(len(coords) - 1):
u = coords[i]
v = coords[i + 1]
self.G.add_edge(u, v, weight=Point(coords[i]).distance(Point(coords[i + 1])))
except Exception as e:
logging.error(f"Error creating graph: {e}")
raise
def _create_mst(self):
"""
Create a Minimum Spanning Tree (MST) from the graph.
"""
def find_paths_between_nearest_points(g: nx.Graph, nearest_points: List[Point]) -> List[Tuple]:
edges = []
for i, start_point in enumerate(nearest_points):
start = (start_point.x, start_point.y)
for end_point in nearest_points[i + 1:]:
end = (end_point.x, end_point.y)
if nx.has_path(g, start, end):
path = nx.shortest_path(g, source=start, target=end, weight='weight')
path_edges = list(zip(path[:-1], path[1:]))
edges.extend((u, v, g[u][v]['weight']) for u, v in path_edges)
return edges
try:
edges = find_paths_between_nearest_points(self.G, self.nearest_points)
h = nx.Graph()
h.add_weighted_edges_from(edges)
mst = nx.minimum_spanning_tree(h, weight='weight')
final_edges = []
for u, v in mst.edges():
if nx.has_path(self.G, u, v):
path = nx.shortest_path(self.G, source=u, target=v, weight='weight')
path_edges = list(zip(path[:-1], path[1:]))
final_edges.extend((x, y, self.G[x][y]['weight']) for x, y in path_edges)
self.final_mst = nx.Graph()
self.final_mst.add_weighted_edges_from(final_edges)
except Exception as e:
logging.error(f"Error creating MST: {e}")
raise
def _iteratively_remove_edges(self):
"""
Iteratively remove edges that do not have any nearest points and have one end with only one connection.
Also remove nodes that don't have any connections and street nodes with only one connection.
"""
nearest_points_tuples = [(point.x, point.y) for point in self.nearest_points]
def find_edges_to_remove(graph: nx.Graph) -> List[Tuple]:
edges_to_remove = []
for u, v, d in graph.edges(data=True):
if u not in nearest_points_tuples and v not in nearest_points_tuples:
if graph.degree(u) == 1 or graph.degree(v) == 1:
edges_to_remove.append((u, v, d))
return edges_to_remove
def find_nodes_to_remove(graph: nx.Graph) -> List[Tuple]:
nodes_to_remove = []
for node in graph.nodes():
if graph.degree(node) == 0:
nodes_to_remove.append(node)
return nodes_to_remove
try:
edges_to_remove = find_edges_to_remove(self.final_mst)
self.final_mst_steps = [list(self.final_mst.edges(data=True))]
while edges_to_remove or find_nodes_to_remove(self.final_mst):
self.final_mst.remove_edges_from(edges_to_remove)
nodes_to_remove = find_nodes_to_remove(self.final_mst)
self.final_mst.remove_nodes_from(nodes_to_remove)
edges_to_remove = find_edges_to_remove(self.final_mst)
self.final_mst_steps.append(list(self.final_mst.edges(data=True)))
def find_single_connection_street_nodes(graph: nx.Graph) -> List[Tuple]:
single_connection_street_nodes = []
for node in graph.nodes():
if node not in nearest_points_tuples and graph.degree(node) == 1:
single_connection_street_nodes.append(node)
return single_connection_street_nodes
single_connection_street_nodes = find_single_connection_street_nodes(self.final_mst)
while single_connection_street_nodes:
for node in single_connection_street_nodes:
neighbors = list(self.final_mst.neighbors(node))
self.final_mst.remove_node(node)
for neighbor in neighbors:
if self.final_mst.degree(neighbor) == 0:
self.final_mst.remove_node(neighbor)
single_connection_street_nodes = find_single_connection_street_nodes(self.final_mst)
self.final_mst_steps.append(list(self.final_mst.edges(data=True)))
except Exception as e:
logging.error(f"Error iteratively removing edges: {e}")
raise
def _add_centroids_to_mst(self):
"""
Add centroids to the final MST graph and connect them to their associated node on the graph.
"""
try:
for i, centroid in enumerate(self.centroids):
building_name = self.building_names[i]
pos = (centroid.x, centroid.y)
node_type = 'building' if 'central_plant' not in building_name else 'generation'
self.final_mst.add_node(pos, type=node_type, name=building_name, pos=pos)
nearest_point = None
min_distance = float('inf')
for node in self.final_mst.nodes():
if self.final_mst.nodes[node].get('type') != 'building' and self.final_mst.nodes[node].get('type') != 'generation':
distance = centroid.distance(Point(node))
if distance < min_distance:
min_distance = distance
nearest_point = node
if nearest_point:
self.final_mst.add_edge(pos, nearest_point, weight=min_distance)
except Exception as e:
logging.error(f"Error adding centroids to MST: {e}")
raise
def _convert_edge_weights_to_meters(self):
"""
Convert all edge weights in the final MST graph to meters using the Haversine formula.
"""
try:
for u, v, data in self.final_mst.edges(data=True):
distance = haversine(u[0], u[1], v[0], v[1])
self.final_mst[u][v]['weight'] = distance
except Exception as e:
logging.error(f"Error converting edge weights to meters: {e}")
raise
def _create_final_network_graph(self):
"""
Create the final network graph with the required attributes from the final MST.
"""
self.network_graph = nx.Graph()
node_id = 1
node_mapping = {}
for node in self.final_mst.nodes:
pos = node
if 'type' in self.final_mst.nodes[node]:
if self.final_mst.nodes[node]['type'] == 'building':
name = self.final_mst.nodes[node]['name']
node_type = 'building'
elif self.final_mst.nodes[node]['type'] == 'generation':
name = self.final_mst.nodes[node]['name']
node_type = 'generation'
else:
name = f'junction_{node_id}'
node_type = 'junction'
self.network_graph.add_node(node_id, name=name, type=node_type, pos=pos)
node_mapping[node] = node_id
node_id += 1
for u, v, data in self.final_mst.edges(data=True):
u_new = node_mapping[u]
v_new = node_mapping[v]
length = data['weight']
self.network_graph.add_edge(u_new, v_new, length=length)
def plot_network_graph(self):
"""
Plot the network graph using matplotlib and networkx.
"""
plt.figure(figsize=(15, 10))
pos = {node: data['pos'] for node, data in self.network_graph.nodes(data=True)}
nx.draw_networkx_nodes(self.network_graph, pos, node_color='blue', node_size=50)
nx.draw_networkx_edges(self.network_graph, pos, edge_color='gray')
plt.title('District Heating Network Graph')
plt.axis('off')
plt.show()

56
scripts/district_heating_network/road_processor.py Normal file
View File

@ -0,0 +1,56 @@
from pathlib import Path
from shapely.geometry import Polygon, Point, shape
import json
def road_processor(x, y, diff):
"""
Processes a .JSON file to find roads that have at least one node within a specified polygon.
Parameters:
x (float): The x-coordinate of the center of the selection box.
y (float): The y-coordinate of the center of the selection box.
diff (float): The half-width of the selection box.
Returns:
str: The file path of the output GeoJSON file containing the selected roads.
"""
diff += 2 * diff
# Define the selection polygon
selection_box = Polygon([
[x + diff, y - diff],
[x - diff, y - diff],
[x - diff, y + diff],
[x + diff, y + diff]
])
# Define input and output file paths
geojson_file = Path("./input_files/roads.json").resolve()
output_file = Path('./input_files/output_roads.geojson').resolve()
# Initialize a list to store the roads in the region
roads_in_region = []
# Read the GeoJSON file
with open(geojson_file, 'r') as file:
roads = json.load(file)
line_features = [feature for feature in roads['features'] if feature['geometry']['type'] == 'LineString']
# Check each road feature
for feature in line_features:
road_shape = shape(feature['geometry'])
# Check if any node of the road is inside the selection box
if any(selection_box.contains(Point(coord)) for coord in road_shape.coords):
roads_in_region.append(feature)
# Create a new GeoJSON structure with the selected roads
output_geojson = {
"type": "FeatureCollection",
"features": roads_in_region
}
# Write the selected roads to the output file
with open(output_file, 'w') as outfile:
json.dump(output_geojson, outfile)
return output_file

54
simulation_result_test.py
View File

@ -1,54 +0,0 @@
from pathlib import Path
import subprocess
from scripts.ep_run_enrich import energy_plus_workflow
from hub.imports.geometry_factory import GeometryFactory
from hub.helpers.dictionaries import Dictionaries
from hub.imports.construction_factory import ConstructionFactory
from hub.imports.usage_factory import UsageFactory
from hub.imports.weather_factory import WeatherFactory
from hub.imports.results_factory import ResultFactory
from scripts.energy_system_retrofit_report import EnergySystemRetrofitReport
from scripts.geojson_creator import process_geojson
from scripts import random_assignation
from hub.imports.energy_systems_factory import EnergySystemsFactory
from scripts.energy_system_sizing import SystemSizing
from scripts.solar_angles import CitySolarAngles
from scripts.pv_sizing_and_simulation import PVSizingSimulation
from scripts.energy_system_retrofit_results import consumption_data, cost_data
from scripts.energy_system_sizing_and_simulation_factory import EnergySystemsSimulationFactory
from scripts.costs.cost import Cost
from scripts.costs.constants import SKIN_RETROFIT_AND_SYSTEM_RETROFIT_AND_PV, SYSTEM_RETROFIT_AND_PV, CURRENT_STATUS
import hub.helpers.constants as cte
from hub.exports.exports_factory import ExportsFactory
from scripts.pv_feasibility import pv_feasibility
# Specify the GeoJSON file path
input_files_path = (Path(__file__).parent / 'input_files')
input_files_path.mkdir(parents=True, exist_ok=True)
geojson_file = process_geojson(x=-73.5681295982132, y=45.49218262677643, diff=0.0001)
geojson_file_path = input_files_path / 'output_buildings.geojson'
output_path = (Path(__file__).parent / 'out_files').resolve()
output_path.mkdir(parents=True, exist_ok=True)
energy_plus_output_path = output_path / 'energy_plus_outputs'
energy_plus_output_path.mkdir(parents=True, exist_ok=True)
simulation_results_path = (Path(__file__).parent / 'out_files' / 'simulation_results').resolve()
simulation_results_path.mkdir(parents=True, exist_ok=True)
sra_output_path = output_path / 'sra_outputs'
sra_output_path.mkdir(parents=True, exist_ok=True)
cost_analysis_output_path = output_path / 'cost_analysis'
cost_analysis_output_path.mkdir(parents=True, exist_ok=True)
city = GeometryFactory(file_type='geojson',
path=geojson_file_path,
height_field='height',
year_of_construction_field='year_of_construction',
function_field='function',
function_to_hub=Dictionaries().montreal_function_to_hub_function).city
ConstructionFactory('nrcan', city).enrich()
UsageFactory('nrcan', city).enrich()
WeatherFactory('epw', city).enrich()
energy_plus_workflow(city, energy_plus_output_path)
random_assignation.call_random(city.buildings, random_assignation.residential_new_systems_percentage)
EnergySystemsFactory('montreal_future', city).enrich()
for building in city.buildings:
EnergySystemsSimulationFactory('archetype1', building=building, output_path=simulation_results_path).enrich()