Codes transfered

.idea/energy_system_modelling_workflow.iml

@@ -2,7 +2,7 @@
 <module type="PYTHON_MODULE" version="4">
   <component name="NewModuleRootManager">
     <content url="file://$MODULE_DIR$" />
-    <orderEntry type="inheritedJdk" />
+    <orderEntry type="jdk" jdkName="hub" jdkType="Python SDK" />
     <orderEntry type="sourceFolder" forTests="false" />
   </component>
 </module>

.idea/misc.xml

@@ -1,4 +1,4 @@
 <?xml version="1.0" encoding="UTF-8"?>
 <project version="4">
-  <component name="ProjectRootManager" version="2" project-jdk-name="Python 3.10" project-jdk-type="Python SDK" />
+  <component name="ProjectRootManager" version="2" project-jdk-name="hub" project-jdk-type="Python SDK" />
 </project>

main.py

@@ -0,0 +1,42 @@
+import os
+
+from scripts.geojson_creator import process_geojson
+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 hub.exports.exports_factory import ExportsFactory
+from scripts.energy_system_analysis_report import EnergySystemAnalysisReport
+from scripts import random_assignation
+from hub.imports.energy_systems_factory import EnergySystemsFactory
+
+# Specify the GeoJSON file path
+geojson_file = process_geojson(x=-73.5681295982132, y=45.49218262677643, diff=0.0001)
+file_path = (Path(__file__).parent.parent / 'input_files' / f'{geojson_file}')
+# Specify the output path for the PDF file
+output_path = (Path(__file__).parent / 'out_files').resolve()
+# Create city object from GeoJSON file
+city = GeometryFactory('geojson',
+                       path=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
+# Enrich city data
+ConstructionFactory('nrcan', city).enrich()
+UsageFactory('comnet', city).enrich()
+WeatherFactory('epw', city).enrich()
+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()
+# Run EnergyPlus workflow
+energy_plus_workflow(city)
+random_assignation.call_random(city.buildings, random_assignation.residential_systems_percentage)
+EnergySystemsFactory('montreal_custom', city).enrich()
+EnergySystemAnalysisReport(city, city.buildings).create_report()

scripts/energy_system_analysis_report.py

@@ -0,0 +1,341 @@
+import subprocess
+import os
+import hub.helpers.constants as cte
+import matplotlib.pyplot as plt
+import random
+import matplotlib.colors as mcolors
+from matplotlib import cm
+import datetime
+from scripts.report_creation import LatexReport
+
+class EnergySystemAnalysisReport:
+  def __init__(self, city, output_path):
+    self.city = city
+    self.output_path = output_path
+    self.content = []
+    self.report = LatexReport('energy_system_analysis_report.tex')
+
+  def building_energy_info(self):
+
+    table_data = [
+      ["Building Name", "Year of Construction", "function", "Yearly Heating Demand (MWh)",
+       "Yearly Cooling Demand (MWh)", "Yearly DHW Demand (MWh)", "Yearly Electricity Demand (MWh)"]
+    ]
+    intensity_table_data = [["Building Name", "Total Floor Area m2", "Heating Demand Intensity kWh/m2",
+                      "Cooling Demand Intensity kWh/m2", "Electricity Intensity kWh/m2"]]
+
+    for building in self.city.buildings:
+      total_floor_area = 0
+      for zone in building.thermal_zones_from_internal_zones:
+        total_floor_area += zone.total_floor_area
+      building_data = [
+        building.name,
+        str(building.year_of_construction),
+        building.function,
+        str(format(building.heating_demand[cte.YEAR][0] / 1e6, '.2f')),
+        str(format(building.cooling_demand[cte.YEAR][0] / 1e6, '.2f')),
+        str(format(building.domestic_hot_water_heat_demand[cte.YEAR][0] / 1e6, '.2f')),
+        str(format((building.lighting_electrical_demand[cte.YEAR][0] + building.appliances_electrical_demand[cte.YEAR][0]) / 1e6, '.2f')),
+      ]
+      intensity_data = [
+        building.name,
+        str(format(total_floor_area, '.2f')),
+        str(format(building.heating_demand[cte.YEAR][0] / (1e3 * total_floor_area), '.2f')),
+        str(format(building.cooling_demand[cte.YEAR][0] / (1e3 * total_floor_area), '.2f')),
+        str(format(
+          (building.lighting_electrical_demand[cte.YEAR][0] + building.appliances_electrical_demand[cte.YEAR][0]) /
+          (1e3 * total_floor_area), '.2f'))
+      ]
+      table_data.append(building_data)
+      intensity_table_data.append(intensity_data)
+
+    self.report.add_table(table_data, caption='City Buildings Energy Demands')
+    self.report.add_table(intensity_table_data, caption='Energy Intensity Information')
+
+  def base_case_charts(self):
+    save_directory = self.output_path
+
+    def create_hvac_demand_chart(building_names, yearly_heating_demand, yearly_cooling_demand):
+      fig, ax = plt.subplots()
+      bar_width = 0.35
+      index = range(len(building_names))
+
+      bars1 = ax.bar(index, yearly_heating_demand, bar_width, label='Yearly Heating Demand (MWh)')
+      bars2 = ax.bar([i + bar_width for i in index], yearly_cooling_demand, bar_width,
+                     label='Yearly Cooling Demand (MWh)')
+
+      ax.set_xlabel('Building Name')
+      ax.set_ylabel('Energy Demand (MWh)')
+      ax.set_title('Yearly HVAC Demands')
+      ax.set_xticks([i + bar_width / 2 for i in index])
+      ax.set_xticklabels(building_names, rotation=45, ha='right')
+      ax.legend()
+      autolabel(bars1, ax)
+      autolabel(bars2, ax)
+      fig.tight_layout()
+      plt.savefig(save_directory / 'hvac_demand_chart.jpg')
+      plt.close()
+
+    def create_bar_chart(title, ylabel, data, filename, bar_color=None):
+      fig, ax = plt.subplots()
+      bar_width = 0.35
+      index = range(len(building_names))
+
+      if bar_color is None:
+        # Generate a random color
+        bar_color = random.choice(list(mcolors.CSS4_COLORS.values()))
+
+      bars = ax.bar(index, data, bar_width, label=ylabel, color=bar_color)
+
+      ax.set_xlabel('Building Name')
+      ax.set_ylabel('Energy Demand (MWh)')
+      ax.set_title(title)
+      ax.set_xticks([i + bar_width / 2 for i in index])
+      ax.set_xticklabels(building_names, rotation=45, ha='right')
+      ax.legend()
+      autolabel(bars, ax)
+      fig.tight_layout()
+      plt.savefig(save_directory / filename)
+      plt.close()
+
+    building_names = [building.name for building in self.city.buildings]
+    yearly_heating_demand = [building.heating_demand[cte.YEAR][0] / 1e6 for building in self.city.buildings]
+    yearly_cooling_demand = [building.cooling_demand[cte.YEAR][0] / 1e6 for building in self.city.buildings]
+    yearly_dhw_demand = [building.domestic_hot_water_heat_demand[cte.YEAR][0] / 1e6 for building in
+                         self.city.buildings]
+    yearly_electricity_demand = [(building.lighting_electrical_demand[cte.YEAR][0] +
+                                  building.appliances_electrical_demand[cte.YEAR][0]) / 1e6 for building in
+                                 self.city.buildings]
+
+    create_hvac_demand_chart(building_names, yearly_heating_demand, yearly_cooling_demand)
+    create_bar_chart('Yearly DHW Demands', 'Energy Demand (MWh)', yearly_dhw_demand, 'dhw_demand_chart.jpg', )
+    create_bar_chart('Yearly Electricity Demands', 'Energy Demand (MWh)', yearly_electricity_demand,
+                     'electricity_demand_chart.jpg')
+    def autolabel(bars, ax):
+      for bar in bars:
+        height = bar.get_height()
+        ax.annotate('{:.1f}'.format(height),
+                    xy=(bar.get_x() + bar.get_width() / 2, height),
+                    xytext=(0, 3),  # 3 points vertical offset
+                    textcoords="offset points",
+                    ha='center', va='bottom')
+
+  def maximum_monthly_hvac_chart(self):
+    save_directory = self.output_path
+    months = ['January', 'February', 'March', 'April', 'May', 'June', 'July', 'August', 'September', 'October',
+              'November', 'December']
+    for building in self.city.buildings:
+      maximum_monthly_heating_load = []
+      maximum_monthly_cooling_load = []
+      fig, axs = plt.subplots(1, 2, figsize=(12, 6))  # Create a figure with 2 subplots side by side
+      for demand in building.heating_peak_load[cte.MONTH]:
+          maximum_monthly_heating_load.append(demand / 3.6e6)
+      for demand in building.cooling_peak_load[cte.MONTH]:
+          maximum_monthly_cooling_load.append(demand / 3.6e6)
+
+    # Plot maximum monthly heating load
+      axs[0].bar(months, maximum_monthly_heating_load, color='red')  # Plot on the first subplot
+      axs[0].set_title('Maximum Monthly Heating Load')
+      axs[0].set_xlabel('Month')
+      axs[0].set_ylabel('Load (kW)')
+      axs[0].tick_params(axis='x', rotation=45)
+
+      # Plot maximum monthly cooling load
+      axs[1].bar(months, maximum_monthly_cooling_load, color='blue')  # Plot on the second subplot
+      axs[1].set_title('Maximum Monthly Cooling Load')
+      axs[1].set_xlabel('Month')
+      axs[1].set_ylabel('Load (kW)')
+      axs[1].tick_params(axis='x', rotation=45)
+
+      plt.tight_layout()  # Adjust layout to prevent overlapping
+      plt.savefig(save_directory / f'{building.name}_monthly_maximum_hvac_loads.jpg')
+      plt.close()
+
+  def load_duration_curves(self):
+    save_directory = self.output_path
+    for building in self.city.buildings:
+      heating_demand = [demand / 1000 for demand in building.heating_demand[cte.HOUR]]
+      cooling_demand = [demand / 1000 for demand in building.cooling_demand[cte.HOUR]]
+      heating_demand_sorted = sorted(heating_demand, reverse=True)
+      cooling_demand_sorted = sorted(cooling_demand, reverse=True)
+
+      plt.style.use('ggplot')
+
+      # Create figure and axis objects with 1 row and 2 columns
+      fig, axs = plt.subplots(1, 2, figsize=(12, 6))
+
+      # Plot sorted heating demand
+      axs[0].plot(heating_demand_sorted, color='red', linewidth=2, label='Heating Demand')
+      axs[0].set_xlabel('Hour', fontsize=14)
+      axs[0].set_ylabel('Heating Demand', fontsize=14)
+      axs[0].set_title('Heating Load Duration Curve', fontsize=16)
+      axs[0].grid(True)
+      axs[0].legend(loc='upper right', fontsize=12)
+
+      # Plot sorted cooling demand
+      axs[1].plot(cooling_demand_sorted, color='blue', linewidth=2, label='Cooling Demand')
+      axs[1].set_xlabel('Hour', fontsize=14)
+      axs[1].set_ylabel('Cooling Demand', fontsize=14)
+      axs[1].set_title('Cooling Load Duration Curve', fontsize=16)
+      axs[1].grid(True)
+      axs[1].legend(loc='upper right', fontsize=12)
+
+      # Adjust layout
+      plt.tight_layout()
+
+      # Save the plot
+      plt.savefig(save_directory / f'{building.name}_load_duration_curve.jpg')
+
+      # Close the plot to release memory
+      plt.close()
+
+  def individual_building_info(self, building):
+    table_data = [
+      ["Maximum Monthly HVAC Demands",
+       f"\\includegraphics[width=1\\linewidth]{{{building.name}_monthly_maximum_hvac_loads.jpg}}"],
+      ["Load Duration Curve", f"\\includegraphics[width=1\\linewidth]{{{building.name}_load_duration_curve.jpg}}"],
+    ]
+
+    self.report.add_table(table_data, caption=f'{building.name} Information', first_column_width=1.5)
+
+  def building_existing_system_info(self, building):
+    existing_archetype = building.energy_systems_archetype_name
+    fuels = []
+    system_schematic = "-"
+    heating_system = "-"
+    cooling_system = "-"
+    dhw = "-"
+    electricity = "Grid"
+    hvac_ec = format((building.heating_consumption[cte.YEAR][0] + building.cooling_consumption[cte.YEAR][0])/1e6, '.2f')
+    dhw_ec = format(building.domestic_hot_water_consumption[cte.YEAR][0]/1e6, '.2f')
+    on_site_generation = "-"
+    yearly_operational_cost = "-"
+    life_cycle_cost = "-"
+    for energy_system in building.energy_systems:
+      if cte.HEATING and cte.DOMESTIC_HOT_WATER in energy_system.demand_types:
+        heating_system = energy_system.name
+        dhw = energy_system.name
+      elif cte.DOMESTIC_HOT_WATER in energy_system.demand_types:
+        dhw = energy_system.name
+      elif cte.HEATING in energy_system.demand_types:
+        heating_system = energy_system.name
+      elif cte.COOLING in energy_system.demand_types:
+        cooling_system = energy_system.name
+      for generation_system in energy_system.generation_systems:
+        fuels.append(generation_system.fuel_type)
+        if generation_system.system_type == cte.PHOTOVOLTAIC:
+          electricity = "Grid-tied PV"
+
+    energy_system_table_data = [
+      ["Detail", "Existing System", "Proposed System"],
+      ["Energy System Archetype", existing_archetype, "-"],
+      ["System Schematic", system_schematic, system_schematic],
+      ["Heating System", heating_system, "-"],
+      ["Cooling System", cooling_system, "-"],
+      ["DHW System", dhw, "-"],
+      ["Electricity", electricity, "-"],
+      ["Fuel(s)", str(fuels), "-"],
+      ["HVAC Energy Consumption (MWh)", hvac_ec, "-"],
+      ["DHW Energy Consumption (MWH)", dhw_ec, "-"],
+      ["Yearly Operational Cost (CAD)", yearly_operational_cost, "-"],
+      ["Life Cycle Cost (CAD)", life_cycle_cost, "-"]
+    ]
+    self.report.add_table(energy_system_table_data, caption= f'Building {building.name} Energy System Characteristics')
+
+  def building_fuel_consumption_breakdown(self, building):
+    save_directory = self.output_path
+    # Initialize variables to store fuel consumption breakdown
+    fuel_breakdown = {
+      "Heating": {"Gas": 0, "Electricity": 0},
+      "Domestic Hot Water": {"Gas": 0, "Electricity": 0},
+      "Cooling": {"Electricity": 0},
+      "Appliance": building.appliances_electrical_demand[cte.YEAR][0] / 1e6,
+      "Lighting": building.lighting_electrical_demand[cte.YEAR][0] / 1e6
+    }
+
+    # Iterate through energy systems of the building
+    for energy_system in building.energy_systems:
+      for demand_type in energy_system.demand_types:
+        for generation_system in energy_system.generation_systems:
+          consumption = 0
+          if demand_type == cte.HEATING:
+            consumption = building.heating_consumption[cte.YEAR][0] / 1e6
+          elif demand_type == cte.DOMESTIC_HOT_WATER:
+            consumption = building.domestic_hot_water_consumption[cte.YEAR][0] / 1e6
+          elif demand_type == cte.COOLING:
+            consumption = building.cooling_consumption[cte.YEAR][0] / 1e6
+
+          if generation_system.fuel_type == cte.ELECTRICITY:
+            fuel_breakdown[demand_type]["Electricity"] += consumption
+          else:
+            fuel_breakdown[demand_type]["Gas"] += consumption
+
+    electricity_labels = ['Appliance', 'Lighting']
+    electricity_sizes = [fuel_breakdown['Appliance'], fuel_breakdown['Lighting']]
+    if fuel_breakdown['Heating']['Electricity'] > 0:
+      electricity_labels.append('Heating')
+      electricity_sizes.append(fuel_breakdown['Heating']['Electricity'])
+    if fuel_breakdown['Cooling']['Electricity'] > 0:
+      electricity_labels.append('Cooling')
+      electricity_sizes.append(fuel_breakdown['Cooling']['Electricity'])
+    if fuel_breakdown['Domestic Hot Water']['Electricity'] > 0:
+      electricity_labels.append('Domestic Hot Water')
+      electricity_sizes.append(fuel_breakdown['Domestic Hot Water']['Electricity'])
+
+    # Data for bar chart
+    gas_labels = ['Heating', 'Domestic Hot Water']
+    gas_sizes = [fuel_breakdown['Heating']['Gas'], fuel_breakdown['Domestic Hot Water']['Gas']]
+
+    # Set the style
+    plt.style.use('ggplot')
+
+    # Create plot grid
+    fig, axs = plt.subplots(1, 2, figsize=(12, 6))
+
+    # Plot pie chart for electricity consumption breakdown
+    colors = cm.get_cmap('tab20c', len(electricity_labels))
+    axs[0].pie(electricity_sizes, labels=electricity_labels,
+               autopct=lambda pct: f"{pct:.1f}%\n({pct / 100 * sum(electricity_sizes):.2f})",
+               startangle=90, colors=[colors(i) for i in range(len(electricity_labels))])
+    axs[0].set_title('Electricity Consumption Breakdown')
+
+    # Plot bar chart for natural gas consumption breakdown
+    colors = cm.get_cmap('Paired', len(gas_labels))
+    axs[1].bar(gas_labels, gas_sizes, color=[colors(i) for i in range(len(gas_labels))])
+    axs[1].set_ylabel('Consumption (MWh)')
+    axs[1].set_title('Natural Gas Consumption Breakdown')
+
+    # Add grid to bar chart
+    axs[1].grid(axis='y', linestyle='--', alpha=0.7)
+
+    # Add a title to the entire figure
+    plt.suptitle('Building Energy Consumption Breakdown', fontsize=16, fontweight='bold')
+
+    # Adjust layout
+    plt.tight_layout()
+
+    # Save the plot as a high-quality image
+    plt.savefig(save_directory / f'{building.name}_energy_consumption_breakdown.png', dpi=300)
+    plt.close()
+
+  def create_report(self):
+    os.chdir(self.output_path)
+    self.report.add_section('Current Status')
+    self.building_energy_info()
+    self.base_case_charts()
+    self.report.add_image('hvac_demand_chart.jpg', caption='Yearly HVAC Demands')
+    self.report.add_image('dhw_demand_chart.jpg', caption='Yearly DHW Demands')
+    self.report.add_image('electricity_demand_chart.jpg', caption='Yearly Electricity Demands')
+    self.maximum_monthly_hvac_chart()
+    self.load_duration_curves()
+    for building in self.city.buildings:
+      self.individual_building_info(building)
+      self.building_existing_system_info(building)
+      self.building_fuel_consumption_breakdown(building)
+      self.report.add_image(f'{building.name}_energy_consumption_breakdown.png',
+                       caption=f'Building {building.name} Consumption by source and sector breakdown')
+    self.report.save_report()
+    self.report.compile_to_pdf()
+
+
+

scripts/energy_system_sizing.py

@@ -0,0 +1,46 @@
+"""
+Energy System rule-based sizing
+SPDX - License - Identifier: LGPL - 3.0 - or -later
+Copyright © 2023 Concordia CERC group
+Project Coder Saeed Ranjbar saeed.ranjbar@concordia.ca
+"""
+
+import hub.helpers.constants as cte
+
+
+class SystemSizing:
+  """
+  The energy system sizing class
+  """
+  def __init__(self, buildings):
+    self.buildings = buildings
+
+  def hvac_sizing(self):
+    for building in self.buildings:
+      peak_heating_demand = building.heating_peak_load[cte.YEAR][0]
+      peak_cooling_demand = building.cooling_peak_load[cte.YEAR][0]
+      if peak_heating_demand > peak_cooling_demand:
+        sizing_demand = peak_heating_demand
+        for system in building.energy_systems:
+          if 'Heating' in system.demand_types:
+            for generation in system.generation_systems:
+              if generation.system_type == 'Heat Pump':
+                generation.nominal_heat_output = 0.6 * sizing_demand / 1000
+                for storage in generation.energy_storage_systems:
+                  if storage.type_energy_stored == 'thermal':
+                    storage.volume = (sizing_demand * 1000) / (cte.WATER_HEAT_CAPACITY*cte.WATER_DENSITY)
+              elif generation.system_type == 'Boiler':
+                generation.nominal_heat_output = 0.4 * sizing_demand / 1000
+
+      else:
+        sizing_demand = peak_cooling_demand
+        for system in building.energy_systems:
+          if 'Cooling' in system.demand_types:
+            for generation in system.generation_systems:
+              if generation.system_type == 'Heat Pump':
+                generation.nominal_heat_output = sizing_demand / 1000
+
+
+
+
+

scripts/ep_run_enrich.py

@@ -0,0 +1,45 @@
+import glob
+import os
+import sys
+from pathlib import Path
+import csv
+from hub.exports.energy_building_exports_factory import EnergyBuildingsExportsFactory
+from hub.imports.results_factory import ResultFactory
+
+sys.path.append('./')
+
+
+def energy_plus_workflow(city):
+  try:
+    # city = city
+    out_path = (Path(__file__).parent.parent / 'out_files')
+    files = glob.glob(f'{out_path}/*')
+
+    for file in files:
+      if file != '.gitignore':
+        os.remove(file)
+    area = 0
+    volume = 0
+    for building in city.buildings:
+      volume = building.volume
+      for ground in building.grounds:
+        area += ground.perimeter_polygon.area
+
+    print('exporting:')
+    _idf = EnergyBuildingsExportsFactory('idf', city, out_path).export()
+    print('    idf exported...')
+    _idf.run()
+
+    csv_file = str((out_path / f'{city.name}_out.csv').resolve())
+    eso_file = str((out_path / f'{city.name}_out.eso').resolve())
+    idf_file = str((out_path / f'{city.name}.idf').resolve())
+    obj_file = str((out_path / f'{city.name}.obj').resolve())
+    ResultFactory('energy_plus_multiple_buildings', city, out_path).enrich()
+
+
+
+  except Exception as ex:
+    print(ex)
+    print('error: ', ex)
+    print('[simulation abort]')
+  sys.stdout.flush()

scripts/geojson_creator.py

@@ -0,0 +1,34 @@
+import json
+from shapely import Polygon
+from shapely import Point
+from pathlib import Path
+
+
+def process_geojson(x, y, diff):
+  selection_box = Polygon([[x + diff, y - diff],
+                           [x - diff, y - diff],
+                           [x - diff, y + diff],
+                           [x + diff, y + diff]])
+  geojson_file = Path('./data/collinear_clean 2.geojson').resolve()
+  output_file = Path('./input_files/output_buildings.geojson').resolve()
+  buildings_in_region = []
+
+  with open(geojson_file, 'r') as file:
+    city = json.load(file)
+    buildings = city['features']
+
+  for building in buildings:
+    coordinates = building['geometry']['coordinates'][0]
+    building_polygon = Polygon(coordinates)
+    centroid = Point(building_polygon.centroid)
+
+    if centroid.within(selection_box):
+      buildings_in_region.append(building)
+
+  output_region = {"type": "FeatureCollection",
+                   "features": buildings_in_region}
+
+  with open(output_file, 'w') as file:
+    file.write(json.dumps(output_region, indent=2))
+
+  return output_file

scripts/random_assignation.py

@@ -0,0 +1,112 @@
+"""
+This project aims to assign energy systems archetype names to Montreal buildings.
+The random assignation is based on statistical information extracted from different sources, being:
+- For residential buildings:
+  - SHEU 2015: https://oee.nrcan.gc.ca/corporate/statistics/neud/dpa/menus/sheu/2015/tables.cfm
+- For non-residential buildings:
+  - Montreal dataportal: https://dataportalforcities.org/north-america/canada/quebec/montreal
+  - https://www.eia.gov/consumption/commercial/data/2018/
+"""
+import json
+import random
+
+from hub.city_model_structure.building import Building
+
+energy_systems_format = 'montreal_custom'
+
+# parameters:
+residential_systems_percentage = {'system 1 gas': 44,
+                                  'system 1 electricity': 6,
+                                  'system 2 gas': 0,
+                                  'system 2 electricity': 0,
+                                  'system 3 and 4 gas': 0,
+                                  'system 3 and 4 electricity': 0,
+                                  'system 5 gas': 0,
+                                  'system 5 electricity': 0,
+                                  'system 6 gas': 0,
+                                  'system 6 electricity': 0,
+                                  'system 8 gas': 44,
+                                  'system 8 electricity': 6}
+
+residential_new_systems_percentage = {'PV+ASHP+GasBoiler+TES': 100,
+                                      'PV+ASHP+ElectricBoiler+TES': 0,
+                                      'PV+GSHP+GasBoiler+TES': 0,
+                                      'PV+GSHP+ElectricBoiler+TES': 0,
+                                      'PV+WSHP+GasBoiler+TES': 0,
+                                      'PV+WSHP+ElectricBoiler+TES': 0}
+
+non_residential_systems_percentage = {'system 1 gas': 0,
+                                      'system 1 electricity': 0,
+                                      'system 2 gas': 0,
+                                      'system 2 electricity': 0,
+                                      'system 3 and 4 gas': 39,
+                                      'system 3 and 4 electricity': 36,
+                                      'system 5 gas': 0,
+                                      'system 5 electricity': 0,
+                                      'system 6 gas': 13,
+                                      'system 6 electricity': 12,
+                                      'system 8 gas': 0,
+                                      'system 8 electricity': 0}
+
+
+def _retrieve_buildings(path, year_of_construction_field=None,
+                        function_field=None, function_to_hub=None, aliases_field=None):
+  _buildings = []
+  with open(path, 'r', encoding='utf8') as json_file:
+    _geojson = json.loads(json_file.read())
+    for feature in _geojson['features']:
+      _building = {}
+      year_of_construction = None
+      if year_of_construction_field is not None:
+        year_of_construction = int(feature['properties'][year_of_construction_field])
+      function = None
+      if function_field is not None:
+        function = feature['properties'][function_field]
+        if function_to_hub is not None:
+          # use the transformation dictionary to retrieve the proper function
+          if function in function_to_hub:
+            function = function_to_hub[function]
+      building_name = ''
+      building_aliases = []
+      if 'id' in feature:
+        building_name = feature['id']
+      if aliases_field is not None:
+        for alias_field in aliases_field:
+          building_aliases.append(feature['properties'][alias_field])
+      _building['year_of_construction'] = year_of_construction
+      _building['function'] = function
+      _building['building_name'] = building_name
+      _building['building_aliases'] = building_aliases
+      _buildings.append(_building)
+  return _buildings
+
+
+def call_random(_buildings: [Building], _systems_percentage):
+  _buildings_with_systems = []
+  _systems_distribution = []
+  _selected_buildings = list(range(0, len(_buildings)))
+  random.shuffle(_selected_buildings)
+  total = 0
+  maximum = 0
+  add_to = 0
+  for _system in _systems_percentage:
+    if _systems_percentage[_system] > 0:
+      number_of_buildings = round(_systems_percentage[_system] / 100 * len(_selected_buildings))
+      _systems_distribution.append({'system': _system, 'number': _systems_percentage[_system],
+                                    'number_of_buildings': number_of_buildings})
+      if number_of_buildings > maximum:
+        maximum = number_of_buildings
+        add_to = len(_systems_distribution) - 1
+      total += number_of_buildings
+  missing = 0
+  if total != len(_selected_buildings):
+    missing = len(_selected_buildings) - total
+  if missing != 0:
+    _systems_distribution[add_to]['number_of_buildings'] += missing
+  _position = 0
+  for case in _systems_distribution:
+    for i in range(0, case['number_of_buildings']):
+      _buildings[_selected_buildings[_position]].energy_systems_archetype_name = case['system']
+      _position += 1
+  return _buildings
+

scripts/report_creation.py

@@ -0,0 +1,78 @@
+import subprocess
+import hub.helpers.constants as cte
+import matplotlib.pyplot as plt
+import random
+import matplotlib.colors as mcolors
+from matplotlib import cm
+import datetime
+
+class LatexReport:
+  def __init__(self, file_name):
+      self.file_name = file_name
+      self.content = []
+      self.content.append(r'\documentclass{article}')
+      self.content.append(r'\usepackage[margin=2.5cm]{geometry}')  # Adjust page margins
+      self.content.append(r'\usepackage{graphicx}')
+      self.content.append(r'\usepackage{tabularx}')
+      self.content.append(r'\begin{document}')
+      # Get current date and time
+      current_datetime = datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S")
+      self.content.append(r'\title{Energy System Analysis Report - ' + current_datetime + r'}')
+      self.content.append(r'\author{Next-Generation Cities Institute}')
+      self.content.append(r'\date{}')  # Remove the date field, as it's included in the title now
+      self.content.append(r'\maketitle')
+
+  def add_section(self, section_title):
+      self.content.append(r'\section{' + section_title + r'}')
+
+  def add_subsection(self, subsection_title):
+      self.content.append(r'\subsection{' + subsection_title + r'}')
+
+  def add_text(self, text):
+      self.content.append(text)
+
+  def add_table(self, table_data, caption=None, first_column_width=None):
+    num_columns = len(table_data[0])
+    total_width = 0.9  # Default total width
+
+    if first_column_width is not None:
+      first_column_width_str = str(first_column_width) + 'cm'
+      total_width -= first_column_width / 16.0  # Adjust total width for the first column
+
+    if caption:
+      self.content.append(r'\begin{table}[htbp]')
+      self.content.append(r'\caption{' + caption + r'}')
+      self.content.append(r'\centering')
+
+    self.content.append(r'\begin{tabularx}{\textwidth}{|p{' + first_column_width_str + r'}|' + '|'.join(['X'] * (
+          num_columns - 1)) + '|}' if first_column_width is not None else r'\begin{tabularx}{\textwidth}{|' + '|'.join(
+      ['X'] * num_columns) + '|}')
+    self.content.append(r'\hline')
+    for row in table_data:
+      self.content.append(' & '.join(row) + r' \\')
+      self.content.append(r'\hline')
+    self.content.append(r'\end{tabularx}')
+
+    if caption:
+      self.content.append(r'\end{table}')
+
+  def add_image(self, image_path, caption=None):
+    if caption:
+      self.content.append(r'\begin{figure}[htbp]')
+      self.content.append(r'\centering')
+      self.content.append(r'\includegraphics[width=0.8\textwidth]{' + image_path + r'}')
+      self.content.append(r'\caption{' + caption + r'}')
+      self.content.append(r'\end{figure}')
+    else:
+      self.content.append(r'\begin{figure}[htbp]')
+      self.content.append(r'\centering')
+      self.content.append(r'\includegraphics[width=0.8\textwidth]{' + image_path + r'}')
+      self.content.append(r'\end{figure}')
+
+  def save_report(self):
+    self.content.append(r'\end{document}')  # Add this line to close the document
+    with open(self.file_name, 'w') as f:
+      f.write('\n'.join(self.content))
+
+  def compile_to_pdf(self):
+      subprocess.run(['pdflatex', self.file_name])

scripts/system_simulation.py

@@ -0,0 +1,132 @@
+import csv
+import math
+from typing import List
+
+from pathlib import Path
+
+import hub.helpers.constants as cte
+from hub.helpers.monthly_values import MonthlyValues
+
+
+class SystemSimulation:
+  def __init__(self, building):
+    self.building = building
+    self.energy_systems = building.energy_systems
+    self.heating_demand = building.heating_demand[cte.HOUR]
+    self.cooling_demand = building.cooling_demand
+    self.dhw_demand = building.domestic_hot_water_heat_demand
+    self.T_out = building.external_temperature[cte.HOUR]
+    self.maximum_heating_demand = building.heating_peak_load[cte.YEAR][0]
+    self.maximum_cooling_demand = building.cooling_peak_load[cte.YEAR][0]
+    self.name = building.name
+    self.energy_system_archetype = building.energy_systems_archetype_name
+
+  def archetype1(self):
+    out_path = (Path(__file__).parent / 'out_files')
+    T, T_sup, T_ret, m_ch, m_dis, q_hp, q_aux = [0] * len(self.heating_demand), [0] * len(
+      self.heating_demand), [0] * len(self.heating_demand), [0] * len(self.heating_demand), [0] * len(
+      self.heating_demand), [0] * len(self.heating_demand), [0] * len(self.heating_demand)
+    hp_electricity: List[float] = [0.0] * len(self.heating_demand)
+    aux_fuel: List[float] = [0.0] * len(self.heating_demand)
+    heating_consumption: List[float] = [0.0] * len(self.heating_demand)
+    boiler_consumption: List[float] = [0.0] * len(self.heating_demand)
+    T[0], dt = 25, 3600  # Assuming dt is defined somewhere
+    ua, v, hp_cap, hp_efficiency, boiler_efficiency = 0, 0, 0, 0, 0
+    for energy_system in self.energy_systems:
+      if cte.ELECTRICITY not in energy_system.demand_types:
+        generation_systems = energy_system.generation_systems
+        for generation_system in generation_systems:
+          if generation_system.system_type == cte.HEAT_PUMP:
+            hp_cap = generation_system.nominal_heat_output
+            hp_efficiency = float(generation_system.heat_efficiency)
+            for storage in generation_system.energy_storage_systems:
+              if storage.type_energy_stored == 'thermal':
+                v, h = float(storage.volume), float(storage.height)
+                r_tot = sum(float(layer.thickness) / float(layer.material.conductivity) for layer in
+                            storage.layers)
+                u_tot = 1 / r_tot
+                d = math.sqrt((4 * v) / (math.pi * h))
+                a_side = math.pi * d * h
+                a_top = math.pi * d ** 2 / 4
+                ua = u_tot * (2 * a_top + a_side)
+          elif generation_system.system_type == cte.BOILER:
+            boiler_cap = generation_system.nominal_heat_output
+            boiler_efficiency = float(generation_system.heat_efficiency)
+
+    for i in range(len(self.heating_demand) - 1):
+      T[i + 1] = T[i] + ((m_ch[i] * (T_sup[i] - T[i])) + (
+              ua * (self.T_out[i] - T[i])) / cte.WATER_HEAT_CAPACITY - m_dis[i] * (T[i] - T_ret[i])) * (dt / (cte.WATER_DENSITY * v))
+      if T[i + 1] < 35:
+          q_hp[i + 1] = hp_cap * 1000
+          m_ch[i + 1] = q_hp[i + 1] / (cte.WATER_HEAT_CAPACITY * 7)
+          T_sup[i + 1] = (q_hp[i + 1] / (m_ch[i + 1] * cte.WATER_HEAT_CAPACITY)) + T[i + 1]
+      elif 35 <= T[i + 1] < 45 and q_hp[i] == 0:
+          q_hp[i + 1] = 0
+          m_ch[i + 1] = 0
+          T_sup[i + 1] = T[i + 1]
+      elif 35 <= T[i + 1] < 45 and q_hp[i] > 0:
+          q_hp[i + 1] = hp_cap * 1000
+          m_ch[i + 1] = q_hp[i + 1] / (cte.WATER_HEAT_CAPACITY * 3)
+          T_sup[i + 1] = (q_hp[i + 1] / (m_ch[i + 1] * cte.WATER_HEAT_CAPACITY)) + T[i + 1]
+      else:
+        q_hp[i + 1], m_ch[i + 1], T_sup[i + 1] = 0, 0, T[i + 1]
+
+      hp_electricity[i + 1] = q_hp[i + 1] / hp_efficiency
+      if self.heating_demand[i + 1] == 0:
+        m_dis[i + 1], t_return, T_ret[i + 1] = 0, T[i + 1], T[i + 1]
+      else:
+        if self.heating_demand[i + 1] > 0.5 * self.maximum_heating_demand:
+          factor = 8
+        else:
+          factor = 4
+        m_dis[i + 1] = self.maximum_heating_demand / (cte.WATER_HEAT_CAPACITY * factor)
+        t_return = T[i + 1] - self.heating_demand[i + 1] / (m_dis[i + 1] * cte.WATER_HEAT_CAPACITY)
+        if m_dis[i + 1] == 0 or (m_dis[i + 1] > 0 and t_return < 25):
+          T_ret[i + 1] = max(25, T[i + 1])
+        else:
+          T_ret[i + 1] = T[i + 1] - self.heating_demand[i + 1] / (m_dis[i + 1] * cte.WATER_HEAT_CAPACITY)
+      tes_output = m_dis[i + 1] * cte.WATER_HEAT_CAPACITY * (T[i + 1] - T_ret[i + 1])
+      if tes_output < self.heating_demand[i + 1]:
+        q_aux[i + 1] = self.heating_demand[i + 1] - tes_output
+        aux_fuel[i + 1] = (q_aux[i + 1] * dt) / 50e6
+      boiler_consumption[i + 1] = q_aux[i + 1] / boiler_efficiency
+      heating_consumption[i + 1] = boiler_consumption[i + 1] + hp_electricity[i + 1]
+    data = list(zip(T, T_sup, T_ret, m_ch, m_dis, q_hp, hp_electricity, aux_fuel, q_aux, self.heating_demand))
+    file_name = f'simulation_results_{self.name}.csv'
+    with open(out_path / file_name, 'w', newline='') as csvfile:
+      output_file = csv.writer(csvfile)
+      # Write header
+      output_file.writerow(['T', 'T_sup', 'T_ret', 'm_ch', 'm_dis', 'q_hp', 'hp_electricity', 'aux_fuel', 'q_aux', 'heating_demand'])
+      # Write data
+      output_file.writerows(data)
+    return heating_consumption, hp_electricity, boiler_consumption
+
+  def enrich(self):
+    if self.energy_system_archetype == 'PV+ASHP+GasBoiler+TES':
+      building_new_heating_consumption, building_heating_electricity_consumption, building_heating_gas_consumption = (
+        self.archetype1())
+      self.building.heating_consumption[cte.HOUR] = building_new_heating_consumption
+      self.building.heating_consumption[cte.MONTH] = MonthlyValues.get_total_month(self.building.heating_consumption[cte.HOUR])
+      self.building.heating_consumption[cte.YEAR] = [sum(self.building.heating_consumption[cte.MONTH])]
+      disaggregated_consumption = {}
+      for energy_system in self.building.energy_systems:
+        if cte.HEATING in energy_system.demand_types:
+          for generation_system in energy_system.generation_systems:
+            disaggregated_consumption[generation_system.fuel_type] = {}
+            if generation_system.fuel_type == cte.ELECTRICITY:
+              disaggregated_consumption[generation_system.fuel_type][
+                cte.HOUR] = building_heating_electricity_consumption
+              disaggregated_consumption[generation_system.fuel_type][cte.MONTH] = MonthlyValues.get_total_month(
+                disaggregated_consumption[generation_system.fuel_type][cte.HOUR])
+              disaggregated_consumption[generation_system.fuel_type][cte.YEAR] = [
+                sum(disaggregated_consumption[generation_system.fuel_type][cte.MONTH])]
+            else:
+              disaggregated_consumption[generation_system.fuel_type][cte.HOUR] = building_heating_gas_consumption
+              disaggregated_consumption[generation_system.fuel_type][cte.MONTH] = MonthlyValues.get_total_month(
+                disaggregated_consumption[generation_system.fuel_type][cte.HOUR])
+              disaggregated_consumption[generation_system.fuel_type][cte.YEAR] = [
+                sum(disaggregated_consumption[generation_system.fuel_type][cte.MONTH])]
+      self.building.heating_consumption_disaggregated = disaggregated_consumption
+    return self.building
+
+