Merge pull request 'OMHM_system_simulation' (#2) from OMHM_system_simulation into main

Reviewed-on: #2

base_case_modelling.py

@@ -0,0 +1,80 @@
+from hub.imports.energy_systems_factory import EnergySystemsFactory
+from hub.imports.results_factory import ResultFactory
+from pathlib import Path
+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
+import json
+import hub.helpers.constants as cte
+from scripts.energy_system_sizing_and_simulation_factory import EnergySystemsSimulationFactory
+
+# Specify the GeoJSON file path
+input_files_path = (Path(__file__).parent / 'input_files')
+geojson_file_path = input_files_path / 'omhm_selected_buildings.geojson'
+output_path = (Path(__file__).parent / 'out_files').resolve()
+output_path.mkdir(parents=True, exist_ok=True)
+ep_output_path = output_path / 'ep_outputs'
+ep_output_path.mkdir(parents=True, exist_ok=True)
+# Create city object from GeoJSON file
+city = GeometryFactory('geojson',
+                       path=geojson_file_path,
+                       height_field='Hieght_LiD',
+                       year_of_construction_field='ANNEE_CONS',
+                       function_field='CODE_UTILI',
+                       function_to_hub=Dictionaries().montreal_function_to_hub_function).city
+# Enrich city data
+ConstructionFactory('nrcan', city).enrich()
+UsageFactory('nrcan', city).enrich()
+WeatherFactory('epw', city).enrich()
+ResultFactory('energy_plus_multiple_buildings', city, ep_output_path).enrich()
+# for building in city.buildings:
+#   building.energy_systems_archetype_name = 'system 7 electricity pv'
+# EnergySystemsFactory('montreal_custom', city).enrich()
+for building in city.buildings:
+  building.energy_systems_archetype_name = 'PV+4Pipe+DHW'
+EnergySystemsFactory('montreal_future', city).enrich()
+for building in city.buildings:
+  EnergySystemsSimulationFactory('archetype13', building=building, output_path=output_path).enrich()
+month_names = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']
+building_data = {}
+for building in city.buildings:
+  building_data[f'building_{building.name}'] = {'id': building.name,
+                                                'total_floor_area':
+                                                  building.thermal_zones_from_internal_zones[0].total_floor_area,
+                                                'yearly_heating_consumption_kWh':
+                                                  building.heating_consumption[cte.YEAR][0] / 3.6e6,
+                                                'yearly_cooling_consumption_kWh':
+                                                  building.cooling_consumption[cte.YEAR][0] / 3.6e6,
+                                                'yearly_dhw_consumption_kWh':
+                                                  building.domestic_hot_water_consumption[cte.YEAR][0] / 3.6e6,
+                                                'yearly_appliance_electricity_consumption_kWh':
+                                                  building.appliances_electrical_demand[cte.YEAR][0] / 3.6e6,
+                                                'yearly_lighting_electricity_consumption_kWh':
+                                                  building.lighting_electrical_demand[cte.YEAR][0] / 3.6e6,
+                                                'heating_peak_load_kW': max(
+                                                  building.heating_consumption[cte.HOUR]) / 3.6e6,
+                                                'cooling_peak_load_kW': max(
+                                                  building.cooling_consumption[cte.HOUR]) / 3.6e6,
+                                                'monthly_heating_demand':
+                                                  {month_name: building.heating_demand[cte.MONTH][i] / 3.6e6
+                                                   for (i, month_name) in enumerate(month_names)},
+                                                'monthly_heating_consumption_kWh':
+                                                  {month_name: building.heating_consumption[cte.MONTH][i] / 3.6e6
+                                                   for (i, month_name) in enumerate(month_names)},
+                                                'monthly_cooling_demand_kWh':
+                                                  {month_name: building.cooling_demand[cte.MONTH][i] / 3.6e6
+                                                   for (i, month_name) in enumerate(month_names)},
+                                                'monthly_cooling_consumption_kWh':
+                                                  {month_name: building.cooling_consumption[cte.MONTH][i] / 3.6e6
+                                                   for (i, month_name) in enumerate(month_names)},
+                                                'monthly_dhw_demand_kWh':
+                                                  {month_name: building.domestic_hot_water_heat_demand[cte.MONTH][i] / 3.6e6
+                                                   for (i, month_name) in enumerate(month_names)},
+                                                'monthly_dhw_consumption_kWh':
+                                                  {month_name: building.domestic_hot_water_consumption[cte.MONTH][i] /
+                                                               3.6e6 for (i, month_name) in enumerate(month_names)}}
+
+with open(output_path / "air_to_water_hp_buildings_data.json", "w") as json_file:
+  json.dump(building_data, json_file, indent=4)

building_modelling/ep_workflow.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('../energy_system_modelling_package/')
+
+
+def energy_plus_workflow(city, output_path):
+  try:
+    # city = city
+    out_path = output_path
+    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()

building_modelling/geojson_creator.py

@@ -0,0 +1,37 @@
+import json
+from shapely import Polygon
+from shapely import Point
+from pathlib import Path
+
+
+def process_geojson(x, y, diff, expansion=False):
+  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()
+  if not expansion:
+    output_file = Path('./input_files/output_buildings.geojson').resolve()
+  else:
+    output_file = Path('./input_files/output_buildings_expanded.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

costing_package/capital_costs.py

@@ -0,0 +1,417 @@
+"""
+Capital costs module
+SPDX - License - Identifier: LGPL - 3.0 - or -later
+Copyright © 2023 Project Coder Guille Gutierrez guillermo.gutierrezmorote@concordia.ca
+Code contributor Pilar Monsalvete Alvarez de Uribarri pilar.monsalvete@concordia.ca
+Code contributor Oriol Gavalda Torrellas oriol.gavalda@concordia.ca
+"""
+import math
+
+import pandas as pd
+import numpy_financial as npf
+from hub.city_model_structure.building import Building
+import hub.helpers.constants as cte
+from costing_package.configuration import Configuration
+from costing_package.constants import (SKIN_RETROFIT, SKIN_RETROFIT_AND_SYSTEM_RETROFIT_AND_PV,
+                             SYSTEM_RETROFIT_AND_PV, CURRENT_STATUS, PV, SYSTEM_RETROFIT)
+from costing_package.cost_base import CostBase
+
+
+class CapitalCosts(CostBase):
+  """
+  Capital costs class
+  """
+
+  def __init__(self, building: Building, configuration: Configuration):
+    super().__init__(building, configuration)
+    self._yearly_capital_costs = pd.DataFrame(
+      index=self._rng,
+      columns=[
+        'B2010_opaque_walls',
+        'B2020_transparent',
+        'B3010_opaque_roof',
+        'B1010_superstructure',
+        'D2010_photovoltaic_system',
+        'D3020_simultaneous_heat_and_cooling_generating_systems',
+        'D3030_heating_systems',
+        'D3040_cooling_systems',
+        'D3050_distribution_systems',
+        'D3060_other_hvac_ahu',
+        'D3070_storage_systems',
+        'D40_dhw',
+      ],
+      dtype='float'
+    )
+    self._yearly_capital_costs.loc[0, 'B2010_opaque_walls'] = 0
+    self._yearly_capital_costs.loc[0, 'B2020_transparent'] = 0
+    self._yearly_capital_costs.loc[0, 'B3010_opaque_roof'] = 0
+    self._yearly_capital_costs.loc[0, 'B1010_superstructure'] = 0
+    self._yearly_capital_costs.loc[0, 'D2010_photovoltaic_system'] = 0
+    self._yearly_capital_costs.loc[0, 'D3020_simultaneous_heat_and_cooling_generating_systems'] = 0
+    self._yearly_capital_costs.loc[0, 'D3030_heating_systems'] = 0
+    self._yearly_capital_costs.loc[0, 'D3040_cooling_systems'] = 0
+    self._yearly_capital_costs.loc[0, 'D3050_distribution_systems'] = 0
+    self._yearly_capital_costs.loc[0, 'D3060_other_hvac_ahu'] = 0
+    self._yearly_capital_costs.loc[0, 'D3070_storage_systems'] = 0
+    self._yearly_capital_costs.loc[0, 'D40_dhw'] = 0
+
+    self._yearly_capital_incomes = pd.DataFrame(
+      index=self._rng,
+      columns=[
+        'Subsidies construction',
+        'Subsidies HVAC',
+        'Subsidies PV'
+      ],
+      dtype='float'
+    )
+    self._yearly_capital_incomes.loc[0, 'Subsidies construction'] = 0
+    self._yearly_capital_incomes.loc[0, 'Subsidies HVAC'] = 0
+    self._yearly_capital_incomes.loc[0, 'Subsidies PV'] = 0
+    self._yearly_capital_costs.fillna(0, inplace=True)
+    self._own_capital = 1 - self._configuration.percentage_credit
+    self._surface_pv = 0
+    for roof in self._building.roofs:
+      self._surface_pv += roof.solid_polygon.area * roof.solar_collectors_area_reduction_factor
+
+    for roof in self._building.roofs:
+      if roof.installed_solar_collector_area is not None:
+        self._surface_pv += roof.installed_solar_collector_area
+      else:
+        self._surface_pv += roof.solid_polygon.area * roof.solar_collectors_area_reduction_factor
+  def calculate(self) -> tuple[pd.DataFrame, pd.DataFrame]:
+    self.skin_capital_cost()
+    self.energy_system_capital_cost()
+    self.skin_yearly_capital_costs()
+    self.yearly_energy_system_costs()
+    self.yearly_incomes()
+    return self._yearly_capital_costs, self._yearly_capital_incomes
+
+  def skin_capital_cost(self):
+    """
+    calculating skin costs
+    :return:
+    """
+    surface_opaque = 0
+    surface_transparent = 0
+    surface_roof = 0
+    surface_ground = 0
+
+    for thermal_zone in self._building.thermal_zones_from_internal_zones:
+      for thermal_boundary in thermal_zone.thermal_boundaries:
+        if thermal_boundary.type == 'Ground':
+          surface_ground += thermal_boundary.opaque_area
+        elif thermal_boundary.type == 'Roof':
+          surface_roof += thermal_boundary.opaque_area
+        elif thermal_boundary.type == 'Wall':
+          surface_opaque += thermal_boundary.opaque_area * (1 - thermal_boundary.window_ratio)
+          surface_transparent += thermal_boundary.opaque_area * thermal_boundary.window_ratio
+
+    chapter = self._capital_costs_chapter.chapter('B_shell')
+    capital_cost_opaque = surface_opaque * chapter.item('B2010_opaque_walls').refurbishment[0]
+    capital_cost_transparent = surface_transparent * chapter.item('B2020_transparent').refurbishment[0]
+    capital_cost_roof = surface_roof * chapter.item('B3010_opaque_roof').refurbishment[0]
+    capital_cost_ground = surface_ground * chapter.item('B1010_superstructure').refurbishment[0]
+    if self._configuration.retrofit_scenario not in (SYSTEM_RETROFIT_AND_PV, CURRENT_STATUS, PV, SYSTEM_RETROFIT):
+      self._yearly_capital_costs.loc[0, 'B2010_opaque_walls'] = capital_cost_opaque * self._own_capital
+      self._yearly_capital_costs.loc[0, 'B2020_transparent'] = capital_cost_transparent * self._own_capital
+      self._yearly_capital_costs.loc[0, 'B3010_opaque_roof'] = capital_cost_roof * self._own_capital
+      self._yearly_capital_costs.loc[0, 'B1010_superstructure'] = capital_cost_ground * self._own_capital
+    capital_cost_skin = capital_cost_opaque + capital_cost_ground + capital_cost_transparent + capital_cost_roof
+    return capital_cost_opaque, capital_cost_transparent, capital_cost_roof, capital_cost_ground, capital_cost_skin
+
+  def skin_yearly_capital_costs(self):
+    skin_capital_cost = self.skin_capital_cost()
+    for year in range(1, self._configuration.number_of_years):
+      self._yearly_capital_costs.loc[year, 'B2010_opaque_walls'] = (
+        -npf.pmt(
+          self._configuration.interest_rate,
+          self._configuration.credit_years,
+          skin_capital_cost[0] * self._configuration.percentage_credit
+        )
+      )
+      self._yearly_capital_costs.loc[year, 'B2020_transparent'] = (
+        -npf.pmt(
+          self._configuration.interest_rate,
+          self._configuration.credit_years,
+          skin_capital_cost[1] * self._configuration.percentage_credit
+        )
+      )
+      self._yearly_capital_costs.loc[year, 'B3010_opaque_roof'] = (
+        -npf.pmt(
+          self._configuration.interest_rate,
+          self._configuration.credit_years,
+          skin_capital_cost[2] * self._configuration.percentage_credit
+        )
+      )
+      self._yearly_capital_costs.loc[year, 'B1010_superstructure'] = (
+        -npf.pmt(
+          self._configuration.interest_rate,
+          self._configuration.credit_years,
+          skin_capital_cost[3] * self._configuration.percentage_credit
+        )
+      )
+
+  def energy_system_capital_cost(self):
+    chapter = self._capital_costs_chapter.chapter('D_services')
+    system_components, component_categories, component_sizes = self.system_components()
+    capital_cost_heating_and_cooling_equipment = 0
+    capital_cost_heating_equipment = 0
+    capital_cost_cooling_equipment = 0
+    capital_cost_domestic_hot_water_equipment = 0
+    capital_cost_energy_storage_equipment = 0
+    capital_cost_distribution_equipment = 0
+    capital_cost_lighting = 0
+    capital_cost_pv = self._surface_pv * chapter.item('D2010_photovoltaic_system').initial_investment[0]
+    for (i, component) in enumerate(system_components):
+      if component_categories[i] == 'multi_generation':
+        capital_cost_heating_and_cooling_equipment += chapter.item(component).initial_investment[0] * component_sizes[i]
+      elif component_categories[i] == 'heating':
+        capital_cost_heating_equipment += chapter.item(component).initial_investment[0] * component_sizes[i]
+      elif component_categories[i] == 'cooling':
+        capital_cost_cooling_equipment += chapter.item(component).initial_investment[0] * component_sizes[i]
+      elif component_categories[i] == 'dhw':
+        capital_cost_domestic_hot_water_equipment += chapter.item(component).initial_investment[0] * \
+                                                      component_sizes[i]
+      elif component_categories[i] == 'distribution':
+        capital_cost_distribution_equipment += chapter.item(component).initial_investment[0] * \
+                                                      component_sizes[i]
+      else:
+        capital_cost_energy_storage_equipment += chapter.item(component).initial_investment[0] * component_sizes[i]
+
+    if self._configuration.retrofit_scenario in (SKIN_RETROFIT_AND_SYSTEM_RETROFIT_AND_PV, SYSTEM_RETROFIT_AND_PV, PV):
+      self._yearly_capital_costs.loc[0, 'D2010_photovoltaic_system'] = capital_cost_pv
+    if (self._configuration.retrofit_scenario in
+       (SYSTEM_RETROFIT_AND_PV, SKIN_RETROFIT_AND_SYSTEM_RETROFIT_AND_PV, SYSTEM_RETROFIT)):
+      self._yearly_capital_costs.loc[0, 'D3020_simultaneous_heat_and_cooling_generating_systems'] = (
+          capital_cost_heating_and_cooling_equipment * self._own_capital)
+      self._yearly_capital_costs.loc[0, 'D3030_heating_systems'] = (
+        capital_cost_heating_equipment * self._own_capital)
+      self._yearly_capital_costs.loc[0, 'D3040_cooling_systems'] = (
+        capital_cost_cooling_equipment * self._own_capital)
+      self._yearly_capital_costs.loc[0, 'D3050_distribution_systems'] = (
+          capital_cost_distribution_equipment * self._own_capital)
+      self._yearly_capital_costs.loc[0, 'D3070_storage_systems'] = (
+          capital_cost_energy_storage_equipment * self._own_capital)
+      self._yearly_capital_costs.loc[0, 'D40_dhw'] = (
+          capital_cost_domestic_hot_water_equipment * self._own_capital)
+    capital_cost_hvac = (capital_cost_heating_and_cooling_equipment + capital_cost_distribution_equipment +
+                         capital_cost_energy_storage_equipment + capital_cost_domestic_hot_water_equipment)
+    return (capital_cost_pv, capital_cost_heating_and_cooling_equipment, capital_cost_heating_equipment,
+            capital_cost_distribution_equipment, capital_cost_cooling_equipment, capital_cost_energy_storage_equipment,
+            capital_cost_domestic_hot_water_equipment, capital_cost_lighting, capital_cost_hvac)
+
+  def yearly_energy_system_costs(self):
+    chapter = self._capital_costs_chapter.chapter('D_services')
+    system_investment_costs = self.energy_system_capital_cost()
+    system_components, component_categories, component_sizes = self.system_components()
+    pv = False
+    for energy_system in self._building.energy_systems:
+      for generation_system in energy_system.generation_systems:
+        if generation_system.system_type == cte.PHOTOVOLTAIC:
+          pv = True
+    for year in range(1, self._configuration.number_of_years):
+      costs_increase = math.pow(1 + self._configuration.consumer_price_index, year)
+      self._yearly_capital_costs.loc[year, 'D2010_photovoltaic_system'] = (
+        -npf.pmt(
+          self._configuration.interest_rate,
+          self._configuration.credit_years,
+          system_investment_costs[0] * self._configuration.percentage_credit
+        )
+      )
+      self._yearly_capital_costs.loc[year, 'D3020_simultaneous_heat_and_cooling_generating_systems'] = (
+        -npf.pmt(
+          self._configuration.interest_rate,
+          self._configuration.credit_years,
+          system_investment_costs[1] * self._configuration.percentage_credit
+        )
+      )
+      self._yearly_capital_costs.loc[year, 'D3030_heating_systems'] = (
+        -npf.pmt(
+          self._configuration.interest_rate,
+          self._configuration.credit_years,
+          system_investment_costs[2] * self._configuration.percentage_credit
+        )
+      )
+      self._yearly_capital_costs.loc[year, 'D3040_cooling_systems'] = (
+        -npf.pmt(
+          self._configuration.interest_rate,
+          self._configuration.credit_years,
+          system_investment_costs[3] * self._configuration.percentage_credit
+        )
+      )
+      self._yearly_capital_costs.loc[year, 'D3050_distribution_systems'] = (
+        -npf.pmt(
+          self._configuration.interest_rate,
+          self._configuration.credit_years,
+          system_investment_costs[4] * self._configuration.percentage_credit
+        )
+      )
+      self._yearly_capital_costs.loc[year, 'D3070_storage_systems'] = (
+        -npf.pmt(
+          self._configuration.interest_rate,
+          self._configuration.credit_years,
+          system_investment_costs[5] * self._configuration.percentage_credit
+        )
+      )
+      self._yearly_capital_costs.loc[year, 'D40_dhw'] = (
+        -npf.pmt(
+          self._configuration.interest_rate,
+          self._configuration.credit_years,
+          system_investment_costs[6] * self._configuration.percentage_credit
+        )
+      )
+      if self._configuration.retrofit_scenario not in (SKIN_RETROFIT, PV):
+        for (i, component) in enumerate(system_components):
+          if (year % chapter.item(component).lifetime) == 0 and year != (self._configuration.number_of_years - 1):
+            if component_categories[i] == 'multi_generation':
+              reposition_cost_heating_and_cooling_equipment = (chapter.item(component).reposition[0] *
+                                                               component_sizes[i] * costs_increase)
+              self._yearly_capital_costs.loc[year, 'D3020_simultaneous_heat_and_cooling_generating_systems'] += (
+                reposition_cost_heating_and_cooling_equipment)
+            elif component_categories[i] == 'heating':
+              reposition_cost_heating_equipment = (chapter.item(component).reposition[0] *
+                                                   component_sizes[i] * costs_increase)
+              self._yearly_capital_costs.loc[year, 'D3030_heating_systems'] += (
+                reposition_cost_heating_equipment)
+            elif component_categories[i] == 'cooling':
+              reposition_cost_cooling_equipment = (chapter.item(component).reposition[0] *
+                                                   component_sizes[i] * costs_increase)
+              self._yearly_capital_costs.loc[year, 'D3040_cooling_systems'] += (
+                reposition_cost_cooling_equipment)
+            elif component_categories[i] == 'dhw':
+              reposition_cost_domestic_hot_water_equipment = (
+                  chapter.item(component).reposition[0] * component_sizes[i] * costs_increase)
+              self._yearly_capital_costs.loc[year, 'D40_dhw'] += reposition_cost_domestic_hot_water_equipment
+            elif component_categories[i] == 'distribution':
+              reposition_cost_distribution_equipment = (
+                  chapter.item(component).reposition[0] * component_sizes[i] * costs_increase)
+              self._yearly_capital_costs.loc[year, 'D3050_distribution_systems'] += (
+                reposition_cost_distribution_equipment)
+            else:
+              reposition_cost_energy_storage_equipment = (
+                  chapter.item(component).initial_investment[0] * component_sizes[i] * costs_increase)
+              self._yearly_capital_costs.loc[year, 'D3070_storage_systems'] += reposition_cost_energy_storage_equipment
+      if self._configuration.retrofit_scenario == CURRENT_STATUS and pv:
+        if (year % chapter.item('D2010_photovoltaic_system').lifetime) == 0:
+          self._yearly_capital_costs.loc[year, 'D2010_photovoltaic_system'] += (
+              self._surface_pv * chapter.item('D2010_photovoltaic_system').reposition[0] * costs_increase
+          )
+      elif self._configuration.retrofit_scenario in (PV, SYSTEM_RETROFIT_AND_PV,
+                                                     SKIN_RETROFIT_AND_SYSTEM_RETROFIT_AND_PV):
+        if (year % chapter.item('D2010_photovoltaic_system').lifetime) == 0:
+          self._yearly_capital_costs.loc[year, 'D2010_photovoltaic_system'] += (
+              self._surface_pv * chapter.item('D2010_photovoltaic_system').reposition[0] * costs_increase
+          )
+
+  def system_components(self):
+    system_components = []
+    component_categories = []
+    sizes = []
+    energy_systems = self._building.energy_systems
+    for energy_system in energy_systems:
+      demand_types = energy_system.demand_types
+      generation_systems = energy_system.generation_systems
+      distribution_systems = energy_system.distribution_systems
+      for generation_system in generation_systems:
+        if generation_system.system_type != cte.PHOTOVOLTAIC:
+          heating_capacity = generation_system.nominal_heat_output or 0
+          cooling_capacity = generation_system.nominal_cooling_output or 0
+          installed_capacity = max(heating_capacity, cooling_capacity) / 1000
+          if cte.DOMESTIC_HOT_WATER in demand_types and cte.HEATING not in demand_types:
+            component_categories.append('dhw')
+            sizes.append(installed_capacity)
+            if generation_system.system_type == cte.HEAT_PUMP:
+              system_components.append(self.heat_pump_type(generation_system, domestic_how_water=True))
+            elif generation_system.system_type == cte.BOILER and generation_system.fuel_type == cte.ELECTRICITY:
+              system_components.append(self.boiler_type(generation_system))
+            else:
+              system_components.append('D302010_template_heat')
+          elif cte.HEATING in demand_types:
+            if cte.COOLING in demand_types and generation_system.fuel_type == cte.ELECTRICITY:
+              component_categories.append('multi_generation')
+            else:
+              component_categories.append('heating')
+            sizes.append(installed_capacity)
+            if generation_system.system_type == cte.HEAT_PUMP:
+              item_type = self.heat_pump_type(generation_system)
+              system_components.append(item_type)
+            elif generation_system.system_type == cte.BOILER:
+              item_type = self.boiler_type(generation_system)
+              system_components.append(item_type)
+            else:
+              if cooling_capacity > heating_capacity:
+                system_components.append('D302090_template_cooling')
+              else:
+                system_components.append('D302010_template_heat')
+          elif cte.COOLING in demand_types:
+            component_categories.append('cooling')
+            sizes.append(installed_capacity)
+            if generation_system.system_type == cte.HEAT_PUMP:
+              item_type = self.heat_pump_type(generation_system)
+              system_components.append(item_type)
+            else:
+              system_components.append('D302090_template_cooling')
+          if generation_system.energy_storage_systems is not None:
+            energy_storage_systems = generation_system.energy_storage_systems
+            for storage_system in energy_storage_systems:
+              if storage_system.type_energy_stored == 'thermal':
+                component_categories.append('thermal storage')
+                sizes.append(storage_system.volume or 0)
+                system_components.append('D306010_storage_tank')
+      if distribution_systems is not None:
+        for distribution_system in distribution_systems:
+          component_categories.append('distribution')
+          sizes.append(self._building.cooling_peak_load[cte.YEAR][0] / 1000)
+          system_components.append('D3040_distribution_systems')
+    return system_components, component_categories, sizes
+
+  def yearly_incomes(self):
+    capital_cost_skin = self.skin_capital_cost()[-1]
+    system_investment_cost = self.energy_system_capital_cost()
+    capital_cost_hvac = system_investment_cost[-1]
+    capital_cost_pv = system_investment_cost[0]
+
+    self._yearly_capital_incomes.loc[0, 'Subsidies construction'] = (
+        capital_cost_skin * self._archetype.income.construction_subsidy/100
+    )
+    self._yearly_capital_incomes.loc[0, 'Subsidies HVAC'] = capital_cost_hvac * self._archetype.income.hvac_subsidy/100
+    self._yearly_capital_incomes.loc[0, 'Subsidies PV'] = capital_cost_pv * self._archetype.income.photovoltaic_subsidy/100
+    self._yearly_capital_incomes.fillna(0, inplace=True)
+
+  @staticmethod
+  def heat_pump_type(generation_system, domestic_how_water=False):
+    source_medium = generation_system.source_medium
+    supply_medium = generation_system.supply_medium
+    if domestic_how_water:
+      heat_pump_item = 'D4010_hot_water_heat_pump'
+    else:
+      if source_medium == cte.AIR and supply_medium == cte.WATER:
+        heat_pump_item = 'D302020_air_to_water_heat_pump'
+      elif source_medium == cte.AIR and supply_medium == cte.AIR:
+        heat_pump_item = 'D302050_air_to_air_heat_pump'
+      elif source_medium == cte.GROUND and supply_medium == cte.WATER:
+        heat_pump_item = 'D302030_ground_to_water_heat_pump'
+      elif source_medium == cte.GROUND and supply_medium == cte.AIR:
+        heat_pump_item = 'D302100_ground_to_air_heat_pump'
+      elif source_medium == cte.WATER and supply_medium == cte.WATER:
+        heat_pump_item = 'D302040_water_to_water_heat_pump'
+      elif source_medium == cte.WATER and supply_medium == cte.AIR:
+        heat_pump_item = 'D302110_water_to_air_heat_pump'
+      else:
+        heat_pump_item = 'D302010_template_heat'
+    return heat_pump_item
+
+  @staticmethod
+  def boiler_type(generation_system):
+    fuel = generation_system.fuel_type
+    if fuel == cte.ELECTRICITY:
+      boiler_item = 'D302080_electrical_boiler'
+    elif fuel == cte.GAS:
+      boiler_item = 'D302070_natural_gas_boiler'
+    else:
+      boiler_item = 'D302010_template_heat'
+    return boiler_item
+
+
+
+

costing_package/configuration.py

@@ -0,0 +1,238 @@
+"""
+Configuration module
+SPDX - License - Identifier: LGPL - 3.0 - or -later
+Copyright © 2023 Project Coder Guille Gutierrez guillermo.gutierrezmorote@concordia.ca
+Code contributor Pilar Monsalvete Alvarez de Uribarri pilar.monsalvete@concordia.ca
+Code contributor Oriol Gavalda Torrellas oriol.gavalda@concordia.ca
+"""
+from hub.catalog_factories.costs_catalog_factory import CostsCatalogFactory
+from hub.catalog_factories.catalog import Catalog
+
+
+class Configuration:
+  """
+  Configuration class
+  """
+
+  def __init__(self,
+               number_of_years,
+               percentage_credit,
+               interest_rate,
+               credit_years,
+               consumer_price_index,
+               electricity_peak_index,
+               electricity_price_index,
+               gas_price_index,
+               discount_rate,
+               retrofitting_year_construction,
+               factories_handler,
+               retrofit_scenario,
+               fuel_type,
+               dictionary,
+               fuel_tariffs
+               ):
+    self._number_of_years = number_of_years
+    self._percentage_credit = percentage_credit
+    self._interest_rate = interest_rate
+    self._credit_years = credit_years
+    self._consumer_price_index = consumer_price_index
+    self._electricity_peak_index = electricity_peak_index
+    self._electricity_price_index = electricity_price_index
+    self._gas_price_index = gas_price_index
+    self._discount_rate = discount_rate
+    self._retrofitting_year_construction = retrofitting_year_construction
+    self._factories_handler = factories_handler
+    self._costs_catalog = CostsCatalogFactory(factories_handler).catalog
+    self._retrofit_scenario = retrofit_scenario
+    self._fuel_type = fuel_type
+    self._dictionary = dictionary
+    self._fuel_tariffs = fuel_tariffs
+
+  @property
+  def number_of_years(self):
+    """
+    Get number of years
+    """
+    return self._number_of_years
+
+  @number_of_years.setter
+  def number_of_years(self, value):
+    """
+    Set number of years
+    """
+    self._number_of_years = value
+
+  @property
+  def percentage_credit(self):
+    """
+    Get percentage credit
+    """
+    return self._percentage_credit
+
+  @percentage_credit.setter
+  def percentage_credit(self, value):
+    """
+    Set percentage credit
+    """
+    self._percentage_credit = value
+
+  @property
+  def interest_rate(self):
+    """
+    Get interest rate
+    """
+    return self._interest_rate
+
+  @interest_rate.setter
+  def interest_rate(self, value):
+    """
+    Set interest rate
+    """
+    self._interest_rate = value
+
+  @property
+  def credit_years(self):
+    """
+    Get credit years
+    """
+    return self._credit_years
+
+  @credit_years.setter
+  def credit_years(self, value):
+    """
+    Set credit years
+    """
+    self._credit_years = value
+
+  @property
+  def consumer_price_index(self):
+    """
+    Get consumer price index
+    """
+    return self._consumer_price_index
+
+  @consumer_price_index.setter
+  def consumer_price_index(self, value):
+    """
+    Set consumer price index
+    """
+    self._consumer_price_index = value
+
+  @property
+  def electricity_peak_index(self):
+    """
+    Get electricity peak index
+    """
+    return self._electricity_peak_index
+
+  @electricity_peak_index.setter
+  def electricity_peak_index(self, value):
+    """
+    Set electricity peak index
+    """
+    self._electricity_peak_index = value
+
+  @property
+  def electricity_price_index(self):
+    """
+    Get electricity price index
+    """
+    return self._electricity_price_index
+
+  @electricity_price_index.setter
+  def electricity_price_index(self, value):
+    """
+    Set electricity price index
+    """
+    self._electricity_price_index = value
+
+  @property
+  def gas_price_index(self):
+    """
+    Get gas price index
+    """
+    return self._gas_price_index
+
+  @gas_price_index.setter
+  def gas_price_index(self, value):
+    """
+    Set gas price index
+    """
+    self._gas_price_index = value
+
+  @property
+  def discount_rate(self):
+    """
+    Get discount rate
+    """
+    return self._discount_rate
+
+  @discount_rate.setter
+  def discount_rate(self, value):
+    """
+    Set discount rate
+    """
+    self._discount_rate = value
+
+  @property
+  def retrofitting_year_construction(self):
+    """
+    Get retrofitting year construction
+    """
+    return self._retrofitting_year_construction
+
+  @retrofitting_year_construction.setter
+  def retrofitting_year_construction(self, value):
+    """
+    Set retrofitting year construction
+    """
+    self._retrofitting_year_construction = value
+
+  @property
+  def factories_handler(self):
+    """
+    Get factories handler
+    """
+    return self._factories_handler
+
+  @factories_handler.setter
+  def factories_handler(self, value):
+    """
+    Set factories handler
+    """
+    self._factories_handler = value
+
+  @property
+  def costs_catalog(self) -> Catalog:
+    """
+    Get costs catalog
+    """
+    return self._costs_catalog
+
+  @property
+  def retrofit_scenario(self):
+    """
+    Get retrofit scenario
+    """
+    return self._retrofit_scenario
+
+  @property
+  def fuel_type(self):
+    """
+    Get fuel type (0: Electricity, 1: Gas)
+    """
+    return self._fuel_type
+
+  @property
+  def dictionary(self):
+    """
+    Get hub function to cost function dictionary
+    """
+    return self._dictionary
+
+  @property
+  def fuel_tariffs(self):
+    """
+    Get fuel tariffs
+    """
+    return self._fuel_tariffs

costing_package/constants.py

@@ -0,0 +1,23 @@
+"""
+Constants module
+SPDX - License - Identifier: LGPL - 3.0 - or -later
+Copyright © 2023 Project Coder Guille Gutierrez guillermo.gutierrezmorote@concordia.ca
+Code contributor Pilar Monsalvete Alvarez de Uribarri pilar.monsalvete@concordia.ca
+Code contributor Oriol Gavalda Torrellas oriol.gavalda@concordia.ca
+"""
+
+# constants
+CURRENT_STATUS = 0
+SKIN_RETROFIT = 1
+SYSTEM_RETROFIT_AND_PV = 2
+SKIN_RETROFIT_AND_SYSTEM_RETROFIT_AND_PV = 3
+PV = 4
+SYSTEM_RETROFIT = 5
+RETROFITTING_SCENARIOS = [
+  CURRENT_STATUS,
+  SKIN_RETROFIT,
+  SYSTEM_RETROFIT_AND_PV,
+  SKIN_RETROFIT_AND_SYSTEM_RETROFIT_AND_PV,
+  PV,
+  SYSTEM_RETROFIT
+]

costing_package/cost.py

@@ -0,0 +1,170 @@
+"""
+Cost module
+SPDX - License - Identifier: LGPL - 3.0 - or -later
+Copyright © 2023 Project Coder Guille Gutierrez guillermo.gutierrezmorote@concordia.ca
+Code contributor Pilar Monsalvete Alvarez de Uribarri pilar.monsalvete@concordia.ca
+Code contributor Oriol Gavalda Torrellas oriol.gavalda@concordia.ca
+"""
+
+import pandas as pd
+import numpy_financial as npf
+from hub.city_model_structure.building import Building
+from hub.helpers.dictionaries import Dictionaries
+from costing_package.configuration import Configuration
+from costing_package.capital_costs import CapitalCosts
+from costing_package.end_of_life_costs import EndOfLifeCosts
+from costing_package.total_maintenance_costs import TotalMaintenanceCosts
+from costing_package.total_operational_costs import TotalOperationalCosts
+from costing_package.total_operational_incomes import TotalOperationalIncomes
+from costing_package.constants import CURRENT_STATUS
+import hub.helpers.constants as cte
+
+
+class Cost:
+  """
+  Cost class
+  """
+
+  def __init__(self,
+               building: Building,
+               number_of_years=31,
+               percentage_credit=0,
+               interest_rate=0.04,
+               credit_years=15,
+               consumer_price_index=0.04,
+               electricity_peak_index=0.05,
+               electricity_price_index=0.05,
+               fuel_price_index=0.05,
+               discount_rate=0.03,
+               retrofitting_year_construction=2020,
+               factories_handler='montreal_new',
+               retrofit_scenario=CURRENT_STATUS,
+               dictionary=None,
+               fuel_tariffs=None):
+    if fuel_tariffs is None:
+      fuel_tariffs = ['Electricity-D', 'Gas-Energir']
+    if dictionary is None:
+      dictionary = Dictionaries().hub_function_to_montreal_custom_costs_function
+    self._building = building
+    fuel_type = self._building.energy_consumption_breakdown.keys()
+    self._configuration = Configuration(number_of_years,
+                                        percentage_credit,
+                                        interest_rate, credit_years,
+                                        consumer_price_index,
+                                        electricity_peak_index,
+                                        electricity_price_index,
+                                        fuel_price_index,
+                                        discount_rate,
+                                        retrofitting_year_construction,
+                                        factories_handler,
+                                        retrofit_scenario,
+                                        fuel_type,
+                                        dictionary,
+                                        fuel_tariffs)
+
+  @property
+  def building(self) -> Building:
+    """
+    Get current building.
+    """
+    return self._building
+
+  def _npv_from_list(self, list_cashflow):
+    return npf.npv(self._configuration.discount_rate, list_cashflow)
+
+  @property
+  def life_cycle(self) -> pd.DataFrame:
+    """
+    Get complete life cycle costs
+    :return: DataFrame
+    """
+    results = pd.DataFrame()
+    global_capital_costs, global_capital_incomes = CapitalCosts(self._building, self._configuration).calculate()
+    global_end_of_life_costs = EndOfLifeCosts(self._building, self._configuration).calculate()
+    global_operational_costs = TotalOperationalCosts(self._building, self._configuration).calculate()
+    global_maintenance_costs = TotalMaintenanceCosts(self._building, self._configuration).calculate()
+    global_operational_incomes = TotalOperationalIncomes(self._building, self._configuration).calculate()
+
+    df_capital_costs_skin = (
+        global_capital_costs['B2010_opaque_walls'] +
+        global_capital_costs['B2020_transparent'] +
+        global_capital_costs['B3010_opaque_roof'] +
+        global_capital_costs['B1010_superstructure']
+    )
+    df_capital_costs_systems = (
+        global_capital_costs['D3020_simultaneous_heat_and_cooling_generating_systems'] +
+        global_capital_costs['D3030_heating_systems'] +
+        global_capital_costs['D3040_cooling_systems'] +
+        global_capital_costs['D3050_distribution_systems'] +
+        global_capital_costs['D3060_other_hvac_ahu'] +
+        global_capital_costs['D3070_storage_systems'] +
+        global_capital_costs['D40_dhw'] +
+        global_capital_costs['D2010_photovoltaic_system']
+    )
+
+    df_end_of_life_costs = global_end_of_life_costs['End_of_life_costs']
+    operational_costs_list = [
+      global_operational_costs['Fixed Costs Electricity Peak'],
+      global_operational_costs['Fixed Costs Electricity Monthly'],
+      global_operational_costs['Variable Costs Electricity']
+      ]
+    additional_costs = [
+                         global_operational_costs[f'Fixed Costs {fuel}'] for fuel in
+                         self._building.energy_consumption_breakdown.keys() if fuel != cte.ELECTRICITY
+                       ] + [
+                         global_operational_costs[f'Variable Costs {fuel}'] for fuel in
+                         self._building.energy_consumption_breakdown.keys() if fuel != cte.ELECTRICITY
+                       ]
+    df_operational_costs = sum(operational_costs_list + additional_costs)
+    df_maintenance_costs = (
+        global_maintenance_costs['Heating_maintenance'] +
+        global_maintenance_costs['Cooling_maintenance'] +
+        global_maintenance_costs['PV_maintenance']
+    )
+    df_operational_incomes = global_operational_incomes['Incomes electricity']
+    df_capital_incomes = (
+        global_capital_incomes['Subsidies construction'] +
+        global_capital_incomes['Subsidies HVAC'] +
+        global_capital_incomes['Subsidies PV']
+    )
+
+    life_cycle_costs_capital_skin = self._npv_from_list(df_capital_costs_skin.values.tolist())
+    life_cycle_costs_capital_systems = self._npv_from_list(df_capital_costs_systems.values.tolist())
+    life_cycle_costs_end_of_life_costs = self._npv_from_list(df_end_of_life_costs.values.tolist())
+    life_cycle_operational_costs = self._npv_from_list(df_operational_costs)
+    life_cycle_maintenance_costs = self._npv_from_list(df_maintenance_costs.values.tolist())
+    life_cycle_operational_incomes = self._npv_from_list(df_operational_incomes.values.tolist())
+    life_cycle_capital_incomes = self._npv_from_list(df_capital_incomes.values.tolist())
+
+    results[f'Scenario {self._configuration.retrofit_scenario}'] = [
+      life_cycle_costs_capital_skin,
+      life_cycle_costs_capital_systems,
+      life_cycle_costs_end_of_life_costs,
+      life_cycle_operational_costs,
+      life_cycle_maintenance_costs,
+      life_cycle_operational_incomes,
+      life_cycle_capital_incomes,
+      global_capital_costs,
+      global_capital_incomes,
+      global_end_of_life_costs,
+      global_operational_costs,
+      global_maintenance_costs,
+      global_operational_incomes
+    ]
+
+    results.index = [
+      'total_capital_costs_skin',
+      'total_capital_costs_systems',
+      'end_of_life_costs',
+      'total_operational_costs',
+      'total_maintenance_costs',
+      'operational_incomes',
+      'capital_incomes',
+      'global_capital_costs',
+      'global_capital_incomes',
+      'global_end_of_life_costs',
+      'global_operational_costs',
+      'global_maintenance_costs',
+      'global_operational_incomes'
+    ]
+    return results

costing_package/cost_base.py

@@ -0,0 +1,40 @@
+"""
+Cost base module
+SPDX - License - Identifier: LGPL - 3.0 - or -later
+Copyright © 2023 Project Coder Guille Gutierrez guillermo.gutierrezmorote@concordia.ca
+Code contributor Pilar Monsalvete Alvarez de Uribarri pilar.monsalvete@concordia.ca
+Code contributor Oriol Gavalda Torrellas oriol.gavalda@concordia.ca
+"""
+
+from hub.city_model_structure.building import Building
+
+from costing_package.configuration import Configuration
+
+
+class CostBase:
+  """
+  Abstract base class for the costs
+  """
+  def __init__(self, building: Building, configuration: Configuration):
+    self._building = building
+    self._configuration = configuration
+    self._total_floor_area = 0
+    for thermal_zone in building.thermal_zones_from_internal_zones:
+      self._total_floor_area += thermal_zone.total_floor_area
+    self._archetype = None
+    self._capital_costs_chapter = None
+    for archetype in self._configuration.costs_catalog.entries().archetypes:
+      if configuration.dictionary[str(building.function)] == str(archetype.function):
+        self._archetype = archetype
+        self._capital_costs_chapter = self._archetype.capital_cost
+        break
+    if not self._archetype:
+      raise KeyError(f'archetype not found for function {building.function}')
+
+    self._rng = range(configuration.number_of_years)
+
+  def calculate(self):
+    """
+    Raises not implemented exception
+    """
+    raise NotImplementedError()

costing_package/end_of_life_costs.py

@@ -0,0 +1,38 @@
+"""
+End of life costs module
+SPDX - License - Identifier: LGPL - 3.0 - or -later
+Copyright © 2023 Project Coder Guille Gutierrez guillermo.gutierrezmorote@concordia.ca
+Code contributor Pilar Monsalvete Alvarez de Uribarri pilar.monsalvete@concordia.ca
+Code contributor Oriol Gavalda Torrellas oriol.gavalda@concordia.ca
+"""
+import math
+import pandas as pd
+from hub.city_model_structure.building import Building
+
+from costing_package.configuration import Configuration
+from costing_package.cost_base import CostBase
+
+
+class EndOfLifeCosts(CostBase):
+  """
+  End of life costs class
+  """
+  def __init__(self, building: Building, configuration: Configuration):
+    super().__init__(building, configuration)
+    self._yearly_end_of_life_costs = pd.DataFrame(index=self._rng, columns=['End_of_life_costs'], dtype='float')
+
+  def calculate(self):
+    """
+    Calculate end of life costs
+    :return: pd.DataFrame
+    """
+    archetype = self._archetype
+    total_floor_area = self._total_floor_area
+    for year in range(1, self._configuration.number_of_years + 1):
+      price_increase = math.pow(1 + self._configuration.consumer_price_index, year)
+      if year == self._configuration.number_of_years:
+        self._yearly_end_of_life_costs.at[year, 'End_of_life_costs'] = (
+            total_floor_area * archetype.end_of_life_cost * price_increase
+        )
+    self._yearly_end_of_life_costs.fillna(0, inplace=True)
+    return self._yearly_end_of_life_costs

costing_package/peak_load.py

@@ -0,0 +1,56 @@
+"""
+Peak load module
+SPDX - License - Identifier: LGPL - 3.0 - or -later
+Copyright © 2023 Project Coder Guille Gutierrez guillermo.gutierrezmorote@concordia.ca
+Code contributor Pilar Monsalvete Alvarez de Uribarri pilar.monsalvete@concordia.ca
+Code contributor Oriol Gavalda Torrellas oriol.gavalda@concordia.ca
+"""
+
+import pandas as pd
+
+import hub.helpers.constants as cte
+
+
+class PeakLoad:
+  """
+  Peak load class
+  """
+
+  def __init__(self, building):
+    self._building = building
+
+  @property
+  def electricity_peak_load(self):
+    """
+    Get the electricity peak load in W
+    """
+    array = [None] * 12
+    heating = 0
+    cooling = 0
+    for system in self._building.energy_systems:
+      if cte.HEATING in system.demand_types:
+        heating = 1
+      if cte.COOLING in system.demand_types:
+        cooling = 1
+    if cte.MONTH in self._building.heating_peak_load.keys() and cte.MONTH in self._building.cooling_peak_load.keys():
+      peak_lighting = self._building.lighting_peak_load[cte.YEAR][0]
+      peak_appliances = self._building.appliances_peak_load[cte.YEAR][0]
+      monthly_electricity_peak = [0.9 * peak_lighting + 0.7 * peak_appliances] * 12
+      conditioning_peak = max(self._building.heating_peak_load[cte.MONTH], self._building.cooling_peak_load[cte.MONTH])
+      for i in range(len(conditioning_peak)):
+        if cooling == 1 and heating == 1:
+          conditioning_peak[i] = conditioning_peak[i]
+          continue
+        elif cooling == 0:
+          conditioning_peak[i] = self._building.heating_peak_load[cte.MONTH][i] * heating
+        else:
+          conditioning_peak[i] = self._building.cooling_peak_load[cte.MONTH][i] * cooling
+        monthly_electricity_peak[i] += 0.8 * conditioning_peak[i]
+
+      electricity_peak_load_results = pd.DataFrame(
+        monthly_electricity_peak,
+        columns=[f'electricity peak load W']
+      )
+    else:
+      electricity_peak_load_results = pd.DataFrame(array, columns=[f'electricity peak load W'])
+    return electricity_peak_load_results

costing_package/total_maintenance_costs.py

@@ -0,0 +1,138 @@
+"""
+Total maintenance costs module
+SPDX - License - Identifier: LGPL - 3.0 - or -later
+Copyright © 2023 Project Coder Guille Gutierrez guillermo.gutierrezmorote@concordia.ca
+Code contributor Pilar Monsalvete Alvarez de Uribarri pilar.monsalvete@concordia.ca
+Code contributor Oriol Gavalda Torrellas oriol.gavalda@concordia.ca
+"""
+import math
+import pandas as pd
+from hub.city_model_structure.building import Building
+import hub.helpers.constants as cte
+
+from costing_package.configuration import Configuration
+from costing_package.cost_base import CostBase
+
+
+class TotalMaintenanceCosts(CostBase):
+  """
+  Total maintenance costs class
+  """
+  def __init__(self, building: Building, configuration: Configuration):
+    super().__init__(building, configuration)
+    self._yearly_maintenance_costs = pd.DataFrame(
+      index=self._rng,
+      columns=[
+        'Heating_maintenance',
+        'Cooling_maintenance',
+        'DHW_maintenance',
+        'PV_maintenance'
+      ],
+      dtype='float'
+    )
+
+  def calculate(self) -> pd.DataFrame:
+    """
+    Calculate total maintenance costs
+    :return: pd.DataFrame
+    """
+    building = self._building
+    archetype = self._archetype
+    # todo: change area pv when the variable exists
+    roof_area = 0
+    surface_pv = 0
+    for roof in self._building.roofs:
+      if roof.installed_solar_collector_area is not None:
+        surface_pv += roof.installed_solar_collector_area
+      else:
+        surface_pv = roof_area * 0.5
+
+    energy_systems = building.energy_systems
+    maintenance_heating_0 = 0
+    maintenance_cooling_0 = 0
+    maintenance_dhw_0 = 0
+    heating_equipments = {}
+    cooling_equipments = {}
+    dhw_equipments = {}
+    for energy_system in energy_systems:
+      if cte.COOLING in energy_system.demand_types:
+        for generation_system in energy_system.generation_systems:
+          if generation_system.fuel_type == cte.ELECTRICITY:
+            if generation_system.system_type == cte.HEAT_PUMP and generation_system.source_medium == cte.AIR:
+              cooling_equipments['air_source_heat_pump'] = generation_system.nominal_cooling_output / 1000
+            elif generation_system.system_type == cte.HEAT_PUMP and generation_system.source_medium == cte.GROUND:
+              cooling_equipments['ground_source_heat_pump'] = generation_system.nominal_cooling_output / 1000
+            elif generation_system.system_type == cte.HEAT_PUMP and generation_system.source_medium == cte.WATER:
+              cooling_equipments['water_source_heat_pump'] = generation_system.nominal_cooling_output / 1000
+            else:
+              cooling_equipments['general_cooling_equipment'] = generation_system.nominal_cooling_output / 1000
+      if cte.HEATING in energy_system.demand_types:
+        for generation_system in energy_system.generation_systems:
+          if generation_system.system_type == cte.HEAT_PUMP and generation_system.source_medium == cte.AIR:
+            heating_equipments['air_source_heat_pump'] = generation_system.nominal_heat_output / 1000
+          elif generation_system.system_type == cte.HEAT_PUMP and generation_system.source_medium == cte.GROUND:
+            heating_equipments['ground_source_heat_pump'] = generation_system.nominal_heat_output / 1000
+          elif generation_system.system_type == cte.HEAT_PUMP and generation_system.source_medium == cte.WATER:
+            heating_equipments['water_source_heat_pump'] = generation_system.nominal_heat_output / 1000
+          elif generation_system.system_type == cte.BOILER and generation_system.fuel_type == cte.GAS:
+            heating_equipments['gas_boiler'] = generation_system.nominal_heat_output / 1000
+          elif generation_system.system_type == cte.BOILER and generation_system.fuel_type == cte.ELECTRICITY:
+            heating_equipments['electric_boiler'] = generation_system.nominal_heat_output / 1000
+          else:
+            heating_equipments['general_heating_equipment'] = generation_system.nominal_heat_output / 1000
+      if cte.DOMESTIC_HOT_WATER in energy_system.demand_types and cte.HEATING not in energy_system.demand_types:
+        for generation_system in energy_system.generation_systems:
+          if generation_system.system_type == cte.HEAT_PUMP and generation_system.source_medium == cte.AIR:
+            dhw_equipments['air_source_heat_pump'] = generation_system.nominal_heat_output / 1000
+          elif generation_system.system_type == cte.HEAT_PUMP and generation_system.source_medium == cte.GROUND:
+            dhw_equipments['ground_source_heat_pump'] = generation_system.nominal_heat_output / 1000
+          elif generation_system.system_type == cte.HEAT_PUMP and generation_system.source_medium == cte.WATER:
+            dhw_equipments['water_source_heat_pump'] = generation_system.nominal_heat_output / 1000
+          elif generation_system.system_type == cte.BOILER and generation_system.fuel_type == cte.GAS:
+            dhw_equipments['gas_boiler'] = generation_system.nominal_heat_output / 1000
+          elif generation_system.system_type == cte.BOILER and generation_system.fuel_type == cte.ELECTRICITY:
+            dhw_equipments['electric_boiler'] = generation_system.nominal_heat_output / 1000
+          else:
+            dhw_equipments['general_heating_equipment'] = generation_system.nominal_heat_output / 1000
+
+
+    for heating_equipment in heating_equipments:
+      component = self.search_hvac_equipment(heating_equipment)
+      maintenance_cost = component.maintenance[0]
+      maintenance_heating_0 += (heating_equipments[heating_equipment] * maintenance_cost)
+
+    for cooling_equipment in cooling_equipments:
+      component = self.search_hvac_equipment(cooling_equipment)
+      maintenance_cost = component.maintenance[0]
+      maintenance_cooling_0 += (cooling_equipments[cooling_equipment] * maintenance_cost)
+
+    for dhw_equipment in dhw_equipments:
+      component = self.search_hvac_equipment(dhw_equipment)
+      maintenance_cost = component.maintenance[0]
+      maintenance_dhw_0 += (dhw_equipments[dhw_equipment] * maintenance_cost)
+
+    maintenance_pv_0 = surface_pv * archetype.operational_cost.maintenance_pv
+
+    for year in range(1, self._configuration.number_of_years + 1):
+      costs_increase = math.pow(1 + self._configuration.consumer_price_index, year)
+      self._yearly_maintenance_costs.loc[year, 'Heating_maintenance'] = (
+          maintenance_heating_0 * costs_increase
+      )
+      self._yearly_maintenance_costs.loc[year, 'Cooling_maintenance'] = (
+          maintenance_cooling_0 * costs_increase
+      )
+      self._yearly_maintenance_costs.loc[year, 'DHW_maintenance'] = (
+          maintenance_dhw_0 * costs_increase
+      )
+      self._yearly_maintenance_costs.loc[year, 'PV_maintenance'] = (
+          maintenance_pv_0 * costs_increase
+      )
+    self._yearly_maintenance_costs.fillna(0, inplace=True)
+    return self._yearly_maintenance_costs
+
+  def search_hvac_equipment(self, equipment_type):
+    for component in self._archetype.operational_cost.maintenance_hvac:
+      if component.type == equipment_type:
+        return component
+
+

costing_package/total_operational_costs.py

@@ -0,0 +1,237 @@
+"""
+Total operational costs module
+SPDX - License - Identifier: LGPL - 3.0 - or -later
+Copyright © 2024 Project Coder Saeed Ranjbar saeed.ranjbar@mail.concordia.ca
+Code contributor Oriol Gavalda Torrellas oriol.gavalda@concordia.ca
+"""
+import math
+import pandas as pd
+
+from hub.city_model_structure.building import Building
+import hub.helpers.constants as cte
+
+from costing_package.configuration import Configuration
+from costing_package.cost_base import CostBase
+from costing_package.peak_load import PeakLoad
+
+
+class TotalOperationalCosts(CostBase):
+  """
+  Total Operational costs class
+  """
+
+  def __init__(self, building: Building, configuration: Configuration):
+    super().__init__(building, configuration)
+    columns_list = self.columns()
+    self._yearly_operational_costs = pd.DataFrame(
+      index=self._rng,
+      columns=columns_list,
+      dtype='float'
+    )
+
+  def calculate(self) -> pd.DataFrame:
+    """
+    Calculate total operational costs
+    :return: pd.DataFrame
+    """
+    building = self._building
+    fuel_consumption_breakdown = building.energy_consumption_breakdown
+    archetype = self._archetype
+    total_floor_area = self._total_floor_area
+    if archetype.function == 'residential':
+      factor = total_floor_area / 80
+    else:
+      factor = 1
+    total_electricity_consumption = sum(self._building.energy_consumption_breakdown[cte.ELECTRICITY].values()) / 3600
+    peak_electricity_load = PeakLoad(self._building).electricity_peak_load
+    peak_load_value = peak_electricity_load.max(axis=1)
+    peak_electricity_demand = peak_load_value[1] / 1000  # self._peak_electricity_demand adapted to kW
+    for system_fuel in self._configuration.fuel_type:
+      fuel = None
+      for fuel_tariff in self._configuration.fuel_tariffs:
+        if system_fuel in fuel_tariff:
+          fuel = self.search_fuel(system_fuel, fuel_tariff)
+      if fuel.type == cte.ELECTRICITY:
+        if fuel.variable.rate_type == 'fixed':
+          variable_electricity_cost_year_0 = (
+              total_electricity_consumption * float(fuel.variable.values[0]) / 1000
+          )
+        else:
+          hourly_electricity_consumption = self.hourly_fuel_consumption_profile(fuel.type)
+          hourly_electricity_price_profile = fuel.variable.values * len(hourly_electricity_consumption)
+          hourly_electricity_price = [hourly_electricity_consumption[i] / 1000 * hourly_electricity_price_profile[i]
+                                      for i in range(len(hourly_electricity_consumption))]
+          variable_electricity_cost_year_0 = sum(hourly_electricity_price)
+        peak_electricity_cost_year_0 = peak_electricity_demand * fuel.fixed_power * 12
+        monthly_electricity_cost_year_0 = fuel.fixed_monthly * 12 * factor
+        for year in range(1, self._configuration.number_of_years + 1):
+          price_increase_electricity = math.pow(1 + self._configuration.electricity_price_index, year)
+          price_increase_peak_electricity = math.pow(1 + self._configuration.electricity_peak_index, year)
+          self._yearly_operational_costs.at[year, 'Fixed Costs Electricity Peak'] = (
+              peak_electricity_cost_year_0 * price_increase_peak_electricity
+          )
+          self._yearly_operational_costs.at[year, 'Fixed Costs Electricity Monthly'] = (
+              monthly_electricity_cost_year_0 * price_increase_peak_electricity
+          )
+          if not isinstance(variable_electricity_cost_year_0, pd.DataFrame):
+            variable_costs_electricity = variable_electricity_cost_year_0 * price_increase_electricity
+          else:
+            variable_costs_electricity = float(variable_electricity_cost_year_0.iloc[0] * price_increase_electricity)
+          self._yearly_operational_costs.at[year, 'Variable Costs Electricity'] = (
+            variable_costs_electricity
+          )
+      else:
+        fuel_fixed_cost = fuel.fixed_monthly * 12 * factor
+        if fuel.type == cte.BIOMASS:
+          conversion_factor = 1
+        else:
+          conversion_factor = fuel.density[0]
+        if fuel.variable.rate_type == 'fixed':
+          variable_cost_fuel = (
+              (sum(fuel_consumption_breakdown[fuel.type].values()) / (
+                    1e6 * fuel.lower_heating_value[0] * conversion_factor)) * fuel.variable.values[0])
+
+        else:
+          hourly_fuel_consumption = self.hourly_fuel_consumption_profile(fuel.type)
+          hourly_fuel_price_profile = fuel.variable.values * len(hourly_fuel_consumption)
+          hourly_fuel_price = [hourly_fuel_consumption[i] / (
+                    1e6 * fuel.lower_heating_value[0] * conversion_factor) * hourly_fuel_price_profile[i]
+                                      for i in range(len(hourly_fuel_consumption))]
+          variable_cost_fuel = sum(hourly_fuel_price)
+
+        for year in range(1, self._configuration.number_of_years + 1):
+          price_increase_gas = math.pow(1 + self._configuration.gas_price_index, year)
+          self._yearly_operational_costs.at[year, f'Fixed Costs {fuel.type}'] = fuel_fixed_cost * price_increase_gas
+          self._yearly_operational_costs.at[year, f'Variable Costs {fuel.type}'] = (
+              variable_cost_fuel * price_increase_gas)
+    self._yearly_operational_costs.fillna(0, inplace=True)
+
+    return self._yearly_operational_costs
+
+  def columns(self):
+    columns_list = []
+    fuels = [key for key in self._building.energy_consumption_breakdown.keys()]
+    for fuel in fuels:
+      if fuel == cte.ELECTRICITY:
+        columns_list.append('Fixed Costs Electricity Peak')
+        columns_list.append('Fixed Costs Electricity Monthly')
+        columns_list.append('Variable Costs Electricity')
+      else:
+        columns_list.append(f'Fixed Costs {fuel}')
+        columns_list.append(f'Variable Costs {fuel}')
+
+    return columns_list
+
+  def search_fuel(self, system_fuel, tariff):
+    fuels = self._archetype.operational_cost.fuels
+    for fuel in fuels:
+      if system_fuel == fuel.type and tariff == fuel.variable.name:
+        return fuel
+    raise KeyError(f'fuel  {system_fuel} with {tariff} tariff not found')
+
+
+  def hourly_fuel_consumption_profile(self, fuel_type):
+    hourly_fuel_consumption = []
+    energy_systems = self._building.energy_systems
+    if fuel_type == cte.ELECTRICITY:
+      appliance = self._building.appliances_electrical_demand[cte.HOUR]
+      lighting = self._building.lighting_electrical_demand[cte.HOUR]
+      elec_heating = 0
+      elec_cooling = 0
+      elec_dhw = 0
+      if cte.HEATING in self._building.energy_consumption_breakdown[cte.ELECTRICITY]:
+        elec_heating = 1
+      if cte.COOLING in self._building.energy_consumption_breakdown[cte.ELECTRICITY]:
+        elec_cooling = 1
+      if cte.DOMESTIC_HOT_WATER in self._building.energy_consumption_breakdown[cte.ELECTRICITY]:
+        elec_dhw = 1
+      heating = None
+      cooling = None
+      dhw = None
+
+      if elec_heating == 1:
+        for energy_system in energy_systems:
+          if cte.HEATING in energy_system.demand_types:
+            for generation_system in energy_system.generation_systems:
+              if generation_system.fuel_type == cte.ELECTRICITY:
+                if cte.HEATING in generation_system.energy_consumption:
+                  heating = generation_system.energy_consumption[cte.HEATING][cte.HOUR]
+                else:
+                  if len(energy_system.generation_systems) > 1:
+                    heating = [x / 2 for x in self._building.heating_consumption[cte.HOUR]]
+                  else:
+                    heating = self._building.heating_consumption[cte.HOUR]
+
+      if elec_dhw == 1:
+        for energy_system in energy_systems:
+          if cte.DOMESTIC_HOT_WATER in energy_system.demand_types:
+            for generation_system in energy_system.generation_systems:
+              if generation_system.fuel_type == cte.ELECTRICITY:
+                if cte.DOMESTIC_HOT_WATER in generation_system.energy_consumption:
+                  dhw = generation_system.energy_consumption[cte.DOMESTIC_HOT_WATER][cte.HOUR]
+                else:
+                  if len(energy_system.generation_systems) > 1:
+                    dhw = [x / 2 for x in self._building.domestic_hot_water_consumption[cte.HOUR]]
+                  else:
+                    dhw = self._building.domestic_hot_water_consumption[cte.HOUR]
+
+      if elec_cooling == 1:
+        for energy_system in energy_systems:
+          if cte.COOLING in energy_system.demand_types:
+            for generation_system in energy_system.generation_systems:
+              if cte.COOLING in generation_system.energy_consumption:
+                cooling = generation_system.energy_consumption[cte.COOLING][cte.HOUR]
+              else:
+                if len(energy_system.generation_systems) > 1:
+                  cooling = [x / 2 for x in self._building.cooling_consumption[cte.HOUR]]
+                else:
+                  cooling = self._building.cooling_consumption[cte.HOUR]
+
+      for i in range(len(self._building.heating_demand[cte.HOUR])):
+        hourly = 0
+        hourly += appliance[i] / 3600
+        hourly += lighting[i] / 3600
+        if heating is not None:
+          hourly += heating[i] / 3600
+        if cooling is not None:
+          hourly += cooling[i] / 3600
+        if dhw is not None:
+          hourly += dhw[i] / 3600
+        hourly_fuel_consumption.append(hourly)
+    else:
+      heating = None
+      dhw = None
+      if cte.HEATING in self._building.energy_consumption_breakdown[fuel_type]:
+        for energy_system in energy_systems:
+          if cte.HEATING in energy_system.demand_types:
+            for generation_system in energy_system.generation_systems:
+              if cte.HEATING in generation_system.energy_consumption:
+                heating = generation_system.energy_consumption[cte.HEATING][cte.HOUR]
+              else:
+                if len(energy_system.generation_systems) > 1:
+                  heating = [x / 2 for x in self._building.heating_consumption[cte.HOUR]]
+                else:
+                  heating = self._building.heating_consumption[cte.HOUR]
+      if cte.DOMESTIC_HOT_WATER in self._building.energy_consumption_breakdown[fuel_type]:
+        for energy_system in energy_systems:
+          if cte.DOMESTIC_HOT_WATER in energy_system.demand_types:
+            for generation_system in energy_system.generation_systems:
+              if cte.DOMESTIC_HOT_WATER in generation_system.energy_consumption:
+                dhw = generation_system.energy_consumption[cte.DOMESTIC_HOT_WATER][cte.HOUR]
+              else:
+                if len(energy_system.generation_systems) > 1:
+                  dhw = [x / 2 for x in self._building.domestic_hot_water_consumption[cte.HOUR]]
+                else:
+                  dhw = self._building.domestic_hot_water_consumption[cte.HOUR]
+
+      for i in range(len(self._building.heating_demand[cte.HOUR])):
+        hourly = 0
+        if heating is not None:
+          hourly += heating[i] / 3600
+        if dhw is not None:
+          hourly += dhw[i] / 3600
+        hourly_fuel_consumption.append(hourly)
+    return hourly_fuel_consumption
+
+
+

costing_package/total_operational_incomes.py

@@ -0,0 +1,44 @@
+"""
+Total operational incomes module
+SPDX - License - Identifier: LGPL - 3.0 - or -later
+Copyright © 2023 Project Coder Guille Gutierrez guillermo.gutierrezmorote@concordia.ca
+Code contributor Pilar Monsalvete Alvarez de Uribarri pilar.monsalvete@concordia.ca
+Code contributor Oriol Gavalda Torrellas oriol.gavalda@concordia.ca
+"""
+import math
+import pandas as pd
+from hub.city_model_structure.building import Building
+import hub.helpers.constants as cte
+
+from costing_package.configuration import Configuration
+from costing_package.cost_base import CostBase
+
+
+class TotalOperationalIncomes(CostBase):
+  """
+  Total operational incomes class
+  """
+  def __init__(self, building: Building, configuration: Configuration):
+    super().__init__(building, configuration)
+    self._yearly_operational_incomes = pd.DataFrame(index=self._rng, columns=['Incomes electricity'], dtype='float')
+
+  def calculate(self) -> pd.DataFrame:
+    """
+    Calculate total operational incomes
+    :return: pd.DataFrame
+    """
+    building = self._building
+    archetype = self._archetype
+    if cte.YEAR not in building.onsite_electrical_production:
+      onsite_electricity_production = 0
+    else:
+      onsite_electricity_production = building.onsite_electrical_production[cte.YEAR][0]
+    for year in range(1, self._configuration.number_of_years + 1):
+      price_increase_electricity = math.pow(1 + self._configuration.electricity_price_index, year)
+      price_export = archetype.income.electricity_export  # to account for unit change
+      self._yearly_operational_incomes.loc[year, 'Incomes electricity'] = (
+        (onsite_electricity_production / 3.6e6) * price_export * price_increase_electricity
+      )
+
+    self._yearly_operational_incomes.fillna(0, inplace=True)
+    return self._yearly_operational_incomes

costing_package/version.py

@@ -0,0 +1,8 @@
+"""
+Cost version number
+SPDX - License - Identifier: LGPL - 3.0 - or -later
+Copyright © 2023 Project Coder Guille Gutierrez guillermo.gutierrezmorote@concordia.ca
+Code contributor Pilar Monsalvete Alvarez de Uribarri pilar.monsalvete@concordia.ca
+Code contributor Oriol Gavalda Torrellas oriol.gavalda@concordia.ca
+"""
+__version__ = '0.1.0.5'

energy_system_modelling_package/energy_system_modelling_factories/archetypes/montreal/archetype_cluster_1.py

@@ -0,0 +1,147 @@
+import csv
+
+from energy_system_modelling_package.energy_system_modelling_factories.hvac_dhw_systems_simulation_models.heat_pump_boiler_tes_heating import \
+  HeatPumpBoilerTesHeating
+from energy_system_modelling_package.energy_system_modelling_factories.hvac_dhw_systems_simulation_models.heat_pump_cooling import \
+  HeatPumpCooling
+from energy_system_modelling_package.energy_system_modelling_factories.hvac_dhw_systems_simulation_models.domestic_hot_water_heat_pump_with_tes import \
+  DomesticHotWaterHeatPumpTes
+from energy_system_modelling_package.energy_system_modelling_factories.pv_assessment.pv_model import PVModel
+from energy_system_modelling_package.energy_system_modelling_factories.pv_assessment.electricity_demand_calculator import HourlyElectricityDemand
+import hub.helpers.constants as cte
+from hub.helpers.monthly_values import MonthlyValues
+
+
+class ArchetypeCluster1:
+  def __init__(self, building, dt, output_path, csv_output=True):
+    self.building = building
+    self.dt = dt
+    self.output_path = output_path
+    self.csv_output = csv_output
+    self.heating_results, self.building_heating_hourly_consumption = self.heating_system_simulation()
+    self.cooling_results, self.total_cooling_consumption_hourly = self.cooling_system_simulation()
+    self.dhw_results, self.total_dhw_consumption_hourly = self.dhw_system_simulation()
+    if 'PV' in self.building.energy_systems_archetype_name:
+      self.pv_results = self.pv_system_simulation()
+    else:
+      self.pv_results = None
+
+  def heating_system_simulation(self):
+    building_heating_hourly_consumption = []
+    boiler = self.building.energy_systems[1].generation_systems[0]
+    hp = self.building.energy_systems[1].generation_systems[1]
+    tes = self.building.energy_systems[1].generation_systems[0].energy_storage_systems[0]
+    heating_demand_joules = self.building.heating_demand[cte.HOUR]
+    heating_peak_load_watts = self.building.heating_peak_load[cte.YEAR][0]
+    upper_limit_tes_heating = 55
+    outdoor_temperature = self.building.external_temperature[cte.HOUR]
+    results = HeatPumpBoilerTesHeating(hp=hp,
+                                       boiler=boiler,
+                                       tes=tes,
+                                       hourly_heating_demand_joules=heating_demand_joules,
+                                       heating_peak_load_watts=heating_peak_load_watts,
+                                       upper_limit_tes=upper_limit_tes_heating,
+                                       outdoor_temperature=outdoor_temperature,
+                                       dt=self.dt).simulation()
+    number_of_ts = int(cte.HOUR_TO_SECONDS/self.dt)
+    heating_consumption_joules = [(x * cte.WATTS_HOUR_TO_JULES) / number_of_ts for x in
+                                  results['Total Heating Power Consumption (W)']]
+    heating_consumption = 0
+    for i in range(1, len(heating_consumption_joules)):
+      heating_consumption += heating_consumption_joules[i]
+      if (i - 1) % number_of_ts == 0:
+        building_heating_hourly_consumption.append(heating_consumption)
+        heating_consumption = 0
+    return results, building_heating_hourly_consumption
+
+  def cooling_system_simulation(self):
+    hp = self.building.energy_systems[2].generation_systems[0]
+    cooling_demand_joules = self.building.cooling_demand[cte.HOUR]
+    cooling_peak_load = self.building.cooling_peak_load[cte.YEAR][0]
+    cutoff_temperature = 11
+    outdoor_temperature = self.building.external_temperature[cte.HOUR]
+    results = HeatPumpCooling(hp=hp,
+                              hourly_cooling_demand_joules=cooling_demand_joules,
+                              cooling_peak_load_watts=cooling_peak_load,
+                              cutoff_temperature=cutoff_temperature,
+                              outdoor_temperature=outdoor_temperature,
+                              dt=self.dt).simulation()
+    building_cooling_hourly_consumption = hp.energy_consumption[cte.COOLING][cte.HOUR]
+    return results, building_cooling_hourly_consumption
+
+  def dhw_system_simulation(self):
+    building_dhw_hourly_consumption = []
+    hp = self.building.energy_systems[-1].generation_systems[0]
+    tes = self.building.energy_systems[-1].generation_systems[0].energy_storage_systems[0]
+    dhw_demand_joules = self.building.domestic_hot_water_heat_demand[cte.HOUR]
+    upper_limit_tes = 65
+    outdoor_temperature = self.building.external_temperature[cte.HOUR]
+    results = DomesticHotWaterHeatPumpTes(hp=hp,
+                                          tes=tes,
+                                          hourly_dhw_demand_joules=dhw_demand_joules,
+                                          upper_limit_tes=upper_limit_tes,
+                                          outdoor_temperature=outdoor_temperature,
+                                          dt=self.dt).simulation()
+    number_of_ts = int(cte.HOUR_TO_SECONDS/self.dt)
+    dhw_consumption_joules = [(x * cte.WATTS_HOUR_TO_JULES) / number_of_ts for x in
+                              results['Total DHW Power Consumption (W)']]
+    dhw_consumption = 0
+    for i in range(1, len(dhw_consumption_joules)):
+      dhw_consumption += dhw_consumption_joules[i]
+      if (i - 1) % number_of_ts == 0:
+        building_dhw_hourly_consumption.append(dhw_consumption)
+        dhw_consumption = 0
+    return results, building_dhw_hourly_consumption
+
+  def pv_system_simulation(self):
+    results = None
+    pv = self.building.energy_systems[0].generation_systems[0]
+    hourly_electricity_demand = HourlyElectricityDemand(self.building).calculate()
+    model_type = 'fixed_efficiency'
+    if model_type == 'fixed_efficiency':
+      results = PVModel(pv=pv,
+                        hourly_electricity_demand_joules=hourly_electricity_demand,
+                        solar_radiation=self.building.roofs[0].global_irradiance_tilted[cte.HOUR],
+                        installed_pv_area=self.building.roofs[0].installed_solar_collector_area,
+                        model_type='fixed_efficiency').fixed_efficiency()
+    return results
+
+  def enrich_building(self):
+    results = self.heating_results | self.cooling_results | self.dhw_results
+    self.building.heating_consumption[cte.HOUR] = self.building_heating_hourly_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])]
+    self.building.cooling_consumption[cte.HOUR] = self.total_cooling_consumption_hourly
+    self.building.cooling_consumption[cte.MONTH] = (
+      MonthlyValues.get_total_month(self.building.cooling_consumption[cte.HOUR]))
+    self.building.cooling_consumption[cte.YEAR] = [sum(self.building.cooling_consumption[cte.MONTH])]
+    self.building.domestic_hot_water_consumption[cte.HOUR] = self.total_dhw_consumption_hourly
+    self.building.domestic_hot_water_consumption[cte.MONTH] = (
+      MonthlyValues.get_total_month(self.building.domestic_hot_water_consumption[cte.HOUR]))
+    self.building.domestic_hot_water_consumption[cte.YEAR] = [
+      sum(self.building.domestic_hot_water_consumption[cte.MONTH])]
+    if self.pv_results is not None:
+      self.building.onsite_electrical_production[cte.HOUR] = [x * cte.WATTS_HOUR_TO_JULES for x in
+                                                              self.pv_results['PV Output (W)']]
+      self.building.onsite_electrical_production[cte.MONTH] = MonthlyValues.get_total_month(self.building.onsite_electrical_production[cte.HOUR])
+      self.building.onsite_electrical_production[cte.YEAR] = [sum(self.building.onsite_electrical_production[cte.MONTH])]
+      if self.csv_output:
+        file_name = f'pv_system_simulation_results_{self.building.name}.csv'
+        with open(self.output_path / file_name, 'w', newline='') as csvfile:
+          output_file = csv.writer(csvfile)
+          # Write header
+          output_file.writerow(self.pv_results.keys())
+          # Write data
+          output_file.writerows(zip(*self.pv_results.values()))
+    if self.csv_output:
+      file_name = f'energy_system_simulation_results_{self.building.name}.csv'
+      with open(self.output_path / file_name, 'w', newline='') as csvfile:
+        output_file = csv.writer(csvfile)
+        # Write header
+        output_file.writerow(results.keys())
+        # Write data
+        output_file.writerows(zip(*results.values()))
+
+
+

energy_system_modelling_package/energy_system_modelling_factories/energy_system_sizing_factory.py

@@ -0,0 +1,79 @@
+"""
+EnergySystemSizingSimulationFactory retrieve the energy system archetype sizing and simulation module
+SPDX - License - Identifier: LGPL - 3.0 - or -later
+Copyright © 2024 Concordia CERC group
+Project Coder Saeed Ranjbar saeed.ranjbar@mail.concordia.ca
+"""
+from energy_system_modelling_package.energy_system_modelling_factories.system_sizing_methods.optimal_sizing import \
+  OptimalSizing
+from energy_system_modelling_package.energy_system_modelling_factories.system_sizing_methods.peak_load_sizing import \
+  PeakLoadSizing
+from energy_system_modelling_package.energy_system_modelling_factories.pv_assessment.pv_sizing import PVSizing
+
+
+class EnergySystemsSizingFactory:
+  """
+  EnergySystemsFactory class
+  """
+
+  def __init__(self, handler, city):
+    self._handler = '_' + handler.lower()
+    self._city = city
+
+  def _peak_load_sizing(self):
+    """
+    Size Energy Systems based on a Load Matching method using the heating, cooling, and dhw peak loads
+    """
+    PeakLoadSizing(self._city).enrich_buildings()
+    self._city.level_of_detail.energy_systems = 1
+    for building in self._city.buildings:
+      building.level_of_detail.energy_systems = 1
+
+  def _optimal_sizing(self):
+    """
+    Size Energy Systems using a Single or Multi Objective GA
+    """
+    OptimalSizing(self._city, optimization_scenario='cost_energy_consumption').enrich_buildings()
+    self._city.level_of_detail.energy_systems = 1
+    for building in self._city.buildings:
+      building.level_of_detail.energy_systems = 1
+
+  def _pv_sizing(self):
+    """
+    Size rooftop, facade or mixture of them for buildings
+    """
+    system_type = 'rooftop'
+    results = {}
+    if system_type == 'rooftop':
+      surface_azimuth = 180
+      maintenance_factor = 0.1
+      mechanical_equipment_factor = 0.3
+      orientation_factor = 0.1
+      tilt_angle = self._city.latitude
+      pv_sizing = PVSizing(self._city,
+                           tilt_angle=tilt_angle,
+                           surface_azimuth=surface_azimuth,
+                           mechanical_equipment_factor=mechanical_equipment_factor,
+                           maintenance_factor=maintenance_factor,
+                           orientation_factor=orientation_factor,
+                           system_type=system_type)
+      results = pv_sizing.rooftop_sizing()
+      pv_sizing.rooftop_tilted_radiation()
+
+    self._city.level_of_detail.energy_systems = 1
+    for building in self._city.buildings:
+      building.level_of_detail.energy_systems = 1
+    return results
+
+  def _district_heating_cooling_sizing(self):
+    """
+    Size District Heating and Cooling Network
+    """
+    pass
+
+  def enrich(self):
+    """
+    Enrich the city given to the class using the class given handler
+    :return: None
+    """
+    return getattr(self, self._handler, lambda: None)()

energy_system_modelling_package/energy_system_modelling_factories/hvac_dhw_systems_simulation_models/domestic_hot_water_heat_pump_with_tes.py

@@ -0,0 +1,120 @@
+import hub.helpers.constants as cte
+from energy_system_modelling_package.energy_system_modelling_factories.hvac_dhw_systems_simulation_models.heat_pump_characteristics import HeatPump
+from energy_system_modelling_package.energy_system_modelling_factories.hvac_dhw_systems_simulation_models.thermal_storage_tank import StorageTank
+from hub.helpers.monthly_values import MonthlyValues
+
+
+class DomesticHotWaterHeatPumpTes:
+  def __init__(self, hp, tes, hourly_dhw_demand_joules, upper_limit_tes,
+               outdoor_temperature, dt=None):
+    self.hp = hp
+    self.tes = tes
+    self.dhw_demand = [demand / cte.WATTS_HOUR_TO_JULES for demand in hourly_dhw_demand_joules]
+    self.upper_limit_tes = upper_limit_tes
+    self.hp_characteristics = HeatPump(self.hp, outdoor_temperature)
+    self.t_out = outdoor_temperature
+    self.dt = dt
+    self.results = {}
+
+  def simulation(self):
+    hp = self.hp
+    tes = self.tes
+    heating_coil_nominal_output = 0
+    if tes.heating_coil_capacity is not None:
+      heating_coil_nominal_output = float(tes.heating_coil_capacity)
+    storage_tank = StorageTank(volume=float(tes.volume),
+                               height=float(tes.height),
+                               material_layers=tes.layers,
+                               heating_coil_capacity=heating_coil_nominal_output)
+
+    hp_delta_t = 8
+    number_of_ts = int(cte.HOUR_TO_SECONDS / self.dt)
+    source_temperature_hourly = self.hp_characteristics.hp_source_temperature()
+    source_temperature = [0] + [x for x in source_temperature_hourly for _ in range(number_of_ts)]
+    demand = [0] + [x for x in self.dhw_demand for _ in range(number_of_ts)]
+    variable_names = ["t_sup_hp", "t_tank", "m_ch", "m_dis", "q_hp", "q_coil", "hp_cop",
+                      "hp_electricity", "available hot water (m3)", "refill flow rate (kg/s)", "total_consumption"]
+    num_hours = len(demand)
+    variables = {name: [0] * num_hours for name in variable_names}
+    (t_sup_hp, t_tank, m_ch, m_dis, m_refill, q_hp, q_coil, hp_cop, hp_electricity, v_dhw, total_consumption) = \
+        [variables[name] for name in variable_names]
+    freshwater_temperature = 18
+    t_tank[0] = 65
+    for i in range(len(demand) - 1):
+      delta_t_demand = demand[i] * (self.dt / (cte.WATER_DENSITY * cte.WATER_HEAT_CAPACITY *
+                                                        storage_tank.volume))
+      if t_tank[i] < self.upper_limit_tes:
+        q_hp[i] = hp.nominal_heat_output
+      delta_t_hp = q_hp[i] * (self.dt / (cte.WATER_DENSITY * cte.WATER_HEAT_CAPACITY * storage_tank.volume))
+      if demand[i] > 0:
+        dhw_needed = (demand[i] * cte.HOUR_TO_SECONDS) / (cte.WATER_HEAT_CAPACITY * t_tank[i] * cte.WATER_DENSITY)
+        m_dis[i] = dhw_needed * cte.WATER_DENSITY / cte.HOUR_TO_SECONDS
+        m_refill[i] = m_dis[i]
+      delta_t_freshwater = m_refill[i] * (t_tank[i] - freshwater_temperature) * (self.dt / (storage_tank.volume *
+                                                                                            cte.WATER_DENSITY))
+      if t_tank[i] < 60:
+        q_coil[i] = float(storage_tank.heating_coil_capacity)
+      delta_t_coil = q_coil[i] * (self.dt / (cte.WATER_DENSITY * cte.WATER_HEAT_CAPACITY * storage_tank.volume))
+      if q_hp[i] > 0:
+        m_ch[i] = q_hp[i] / (cte.WATER_HEAT_CAPACITY * hp_delta_t)
+        t_sup_hp[i] = (q_hp[i] / (m_ch[i] * cte.WATER_HEAT_CAPACITY)) + t_tank[i]
+      else:
+        m_ch[i] = 0
+        t_sup_hp[i] = t_tank[i]
+      if q_hp[i] > 0:
+        if hp.source_medium == cte.AIR and hp.supply_medium == cte.WATER:
+          hp_cop[i] = self.hp_characteristics.air_to_water_cop(source_temperature[i], t_tank[i],
+                                                               mode=cte.DOMESTIC_HOT_WATER)
+        hp_electricity[i] = q_hp[i] / hp_cop[i]
+      else:
+        hp_cop[i] = 0
+        hp_electricity[i] = 0
+
+      t_tank[i + 1] = t_tank[i] + (delta_t_hp - delta_t_freshwater - delta_t_demand + delta_t_coil)
+      total_consumption[i] = hp_electricity[i] + q_coil[i]
+    tes.temperature = []
+    hp_electricity_j = [(x * cte.WATTS_HOUR_TO_JULES) / number_of_ts for x in hp_electricity]
+    heating_coil_j = [(x * cte.WATTS_HOUR_TO_JULES) / number_of_ts for x in q_coil]
+    hp_hourly = []
+    coil_hourly = []
+    coil_sum = 0
+    hp_sum = 0
+    for i in range(1, len(demand)):
+      hp_sum += hp_electricity_j[i]
+      coil_sum += heating_coil_j[i]
+      if (i - 1) % number_of_ts == 0:
+        tes.temperature.append(t_tank[i])
+        hp_hourly.append(hp_sum)
+        coil_hourly.append(coil_sum)
+        hp_sum = 0
+        coil_sum = 0
+    hp.energy_consumption[cte.DOMESTIC_HOT_WATER] = {}
+    hp.energy_consumption[cte.DOMESTIC_HOT_WATER][cte.HOUR] = hp_hourly
+    if len(self.dhw_demand) == 8760:
+      hp.energy_consumption[cte.DOMESTIC_HOT_WATER][cte.MONTH] = MonthlyValues.get_total_month(
+        hp.energy_consumption[cte.DOMESTIC_HOT_WATER][cte.HOUR])
+      hp.energy_consumption[cte.DOMESTIC_HOT_WATER][cte.YEAR] = [
+        sum(hp.energy_consumption[cte.DOMESTIC_HOT_WATER][cte.MONTH])]
+    if self.tes.heating_coil_capacity is not None:
+      tes.heating_coil_energy_consumption[cte.DOMESTIC_HOT_WATER] = {}
+      tes.heating_coil_energy_consumption[cte.DOMESTIC_HOT_WATER][cte.HOUR] = coil_hourly
+      if len(self.dhw_demand) == 8760:
+        tes.heating_coil_energy_consumption[cte.DOMESTIC_HOT_WATER][cte.MONTH] = MonthlyValues.get_total_month(
+          tes.heating_coil_energy_consumption[cte.DOMESTIC_HOT_WATER][cte.HOUR])
+        tes.heating_coil_energy_consumption[cte.DOMESTIC_HOT_WATER][cte.YEAR] = [
+          sum(tes.heating_coil_energy_consumption[cte.DOMESTIC_HOT_WATER][cte.MONTH])]
+
+    self.results['DHW Demand (W)'] = demand
+    self.results['DHW HP Heat Output (W)'] = q_hp
+    self.results['DHW HP Electricity Consumption (W)'] = hp_electricity
+    self.results['DHW HP Source Temperature'] = source_temperature
+    self.results['DHW HP Supply Temperature'] = t_sup_hp
+    self.results['DHW HP COP'] = hp_cop
+    self.results['DHW TES Heating Coil Heat Output (W)'] = q_coil
+    self.results['DHW TES Temperature'] = t_tank
+    self.results['DHW TES Charging Flow Rate (kg/s)'] = m_ch
+    self.results['DHW Flow Rate (kg/s)'] = m_dis
+    self.results['DHW TES Refill Flow Rate (kg/s)'] = m_refill
+    self.results['Available Water in Tank (m3)'] = v_dhw
+    self.results['Total DHW Power Consumption (W)'] = total_consumption
+    return self.results

energy_system_modelling_package/energy_system_modelling_factories/hvac_dhw_systems_simulation_models/heat_pump_boiler_tes_heating.py

@@ -0,0 +1,171 @@
+import hub.helpers.constants as cte
+from hub.helpers.monthly_values import MonthlyValues
+from energy_system_modelling_package.energy_system_modelling_factories.hvac_dhw_systems_simulation_models.heat_pump_characteristics import \
+  HeatPump
+from energy_system_modelling_package.energy_system_modelling_factories.hvac_dhw_systems_simulation_models.thermal_storage_tank import \
+  StorageTank
+
+
+class HeatPumpBoilerTesHeating:
+  def __init__(self, hp, boiler, tes, hourly_heating_demand_joules, heating_peak_load_watts, upper_limit_tes,
+               outdoor_temperature, dt=None):
+    self.hp = hp
+    self.boiler = boiler
+    self.tes = tes
+    self.heating_demand = [demand / cte.WATTS_HOUR_TO_JULES for demand in hourly_heating_demand_joules]
+    if heating_peak_load_watts is not None:
+      self.heating_peak_load = heating_peak_load_watts
+    else:
+      self.heating_peak_load = max(hourly_heating_demand_joules) / cte.HOUR_TO_SECONDS
+    self.upper_limit_tes = upper_limit_tes
+    self.hp_characteristics = HeatPump(self.hp, outdoor_temperature)
+    self.t_out = outdoor_temperature
+    self.dt = dt
+    self.results = {}
+
+  def simulation(self):
+    hp, boiler, tes = self.hp, self.boiler, self.tes
+    heating_coil_nominal_output = 0
+    if tes.heating_coil_capacity is not None:
+      heating_coil_nominal_output = float(tes.heating_coil_capacity)
+    storage_tank = StorageTank(volume=float(tes.volume),
+                               height=float(tes.height),
+                               material_layers=tes.layers,
+                               heating_coil_capacity=heating_coil_nominal_output)
+    number_of_ts = int(cte.HOUR_TO_SECONDS / self.dt)
+    demand = [0] + [x for x in self.heating_demand for _ in range(number_of_ts)]
+    t_out = [0] + [x for x in self.t_out for _ in range(number_of_ts)]
+    source_temperature_hourly = self.hp_characteristics.hp_source_temperature()
+    source_temperature = [0] + [x for x in source_temperature_hourly for _ in range(number_of_ts)]
+    variable_names = ["t_sup_hp", "t_tank", "t_ret", "m_ch", "m_dis", "q_hp", "q_boiler", "hp_cop",
+                      "hp_electricity", "boiler_gas_consumption", "t_sup_boiler", "boiler_energy_consumption",
+                      "heating_consumption", "heating_coil_output", "total_heating_energy_consumption"]
+    num_hours = len(demand)
+    variables = {name: [0] * num_hours for name in variable_names}
+    (t_sup_hp, t_tank, t_ret, m_ch, m_dis, q_hp, q_boiler, hp_cop,
+     hp_electricity, boiler_fuel_consumption, t_sup_boiler, boiler_energy_consumption, heating_consumption, q_coil,
+     total_consumption) = [variables[name] for name in variable_names]
+    t_tank[0] = self.upper_limit_tes
+    hp_delta_t = 5
+    # storage temperature prediction
+    for i in range(len(demand) - 1):
+      t_tank[i + 1] = storage_tank.calculate_space_heating_fully_mixed(charging_flow_rate=m_ch[i],
+                                                                       discharging_flow_rate=m_dis[i],
+                                                                       supply_temperature=t_sup_boiler[i],
+                                                                       return_temperature=t_ret[i],
+                                                                       current_tank_temperature=t_tank[i],
+                                                                       heat_generator_input=q_coil[i],
+                                                                       ambient_temperature=t_out[i],
+                                                                       dt=self.dt)
+      # hp operation
+      if t_tank[i + 1] < 40:
+        q_hp[i + 1] = hp.nominal_heat_output
+        m_ch[i + 1] = q_hp[i + 1] / (cte.WATER_HEAT_CAPACITY * hp_delta_t)
+        t_sup_hp[i + 1] = (q_hp[i + 1] / (m_ch[i + 1] * cte.WATER_HEAT_CAPACITY)) + t_tank[i + 1]
+      elif 40 <= t_tank[i + 1] < self.upper_limit_tes and q_hp[i] == 0:
+        q_hp[i + 1] = 0
+        m_ch[i + 1] = 0
+        t_sup_hp[i + 1] = t_tank[i + 1]
+      elif 40 <= t_tank[i + 1] < self.upper_limit_tes and q_hp[i] > 0:
+        q_hp[i + 1] = hp.nominal_heat_output
+        m_ch[i + 1] = q_hp[i + 1] / (cte.WATER_HEAT_CAPACITY * hp_delta_t)
+        t_sup_hp[i + 1] = (q_hp[i + 1] / (m_ch[i + 1] * cte.WATER_HEAT_CAPACITY)) + t_tank[i + 1]
+      else:
+        q_hp[i + 1], m_ch[i + 1], t_sup_hp[i + 1] = 0, 0, t_tank[i + 1]
+      if q_hp[i + 1] > 0:
+        if hp.source_medium == cte.AIR and self.hp.supply_medium == cte.WATER:
+          hp_cop[i + 1] = self.hp_characteristics.air_to_water_cop(source_temperature[i + 1], t_tank[i + 1], mode=cte.HEATING)
+          hp_electricity[i + 1] = q_hp[i + 1] / hp_cop[i + 1]
+      else:
+        hp_cop[i + 1] = 0
+        hp_electricity[i + 1] = 0
+      # boiler operation
+      if q_hp[i + 1] > 0:
+        if t_sup_hp[i + 1] < 45:
+          q_boiler[i + 1] = boiler.nominal_heat_output
+        elif demand[i + 1] > 0.5 * self.heating_peak_load / self.dt:
+          q_boiler[i + 1] = 0.5 * boiler.nominal_heat_output
+        boiler_energy_consumption[i + 1] = q_boiler[i + 1] / float(boiler.heat_efficiency)
+        if boiler.fuel_type == cte.ELECTRICITY:
+          boiler_fuel_consumption[i + 1] = boiler_energy_consumption[i + 1]
+        else:
+            # TODO: Other fuels should be considered
+          boiler_fuel_consumption[i + 1] = (q_boiler[i + 1] * self.dt) / (
+              float(boiler.heat_efficiency) * cte.NATURAL_GAS_LHV)
+        t_sup_boiler[i + 1] = t_sup_hp[i + 1] + (q_boiler[i + 1] / (m_ch[i + 1] * cte.WATER_HEAT_CAPACITY))
+      # heating coil operation
+      if t_tank[i + 1] < 35:
+        q_coil[i + 1] = heating_coil_nominal_output
+      # storage discharging
+      if demand[i + 1] == 0:
+        m_dis[i + 1] = 0
+        t_ret[i + 1] = t_tank[i + 1]
+      else:
+        if demand[i + 1] > 0.5 * self.heating_peak_load:
+          factor = 8
+        else:
+          factor = 4
+        m_dis[i + 1] = self.heating_peak_load / (cte.WATER_HEAT_CAPACITY * factor)
+        t_ret[i + 1] = t_tank[i + 1] - demand[i + 1] / (m_dis[i + 1] * cte.WATER_HEAT_CAPACITY)
+      # total consumption
+      total_consumption[i + 1] = hp_electricity[i + 1] + boiler_energy_consumption[i + 1] + q_coil[i + 1]
+    tes.temperature = []
+    hp_electricity_j = [(x * cte.WATTS_HOUR_TO_JULES) / number_of_ts for x in hp_electricity]
+    boiler_consumption_j = [(x * cte.WATTS_HOUR_TO_JULES) / number_of_ts for x in boiler_energy_consumption]
+    heating_coil_j = [(x * cte.WATTS_HOUR_TO_JULES) / number_of_ts for x in q_coil]
+    hp_hourly = []
+    boiler_hourly = []
+    coil_hourly = []
+    boiler_sum = 0
+    hp_sum = 0
+    coil_sum = 0
+    for i in range(1, len(demand)):
+      hp_sum += hp_electricity_j[i]
+      boiler_sum += boiler_consumption_j[i]
+      coil_sum += heating_coil_j[i]
+      if (i - 1) % number_of_ts == 0:
+        tes.temperature.append(t_tank[i])
+        hp_hourly.append(hp_sum)
+        boiler_hourly.append(boiler_sum)
+        coil_hourly.append(coil_sum)
+        hp_sum = 0
+        boiler_sum = 0
+        coil_sum = 0
+    hp.energy_consumption[cte.HEATING] = {}
+    hp.energy_consumption[cte.HEATING][cte.HOUR] = hp_hourly
+    boiler.energy_consumption[cte.HEATING] = {}
+    boiler.energy_consumption[cte.HEATING][cte.HOUR] = boiler_hourly
+    if len(self.heating_demand) == 8760:
+      hp.energy_consumption[cte.HEATING][cte.MONTH] = MonthlyValues.get_total_month(
+        hp.energy_consumption[cte.HEATING][cte.HOUR])
+      hp.energy_consumption[cte.HEATING][cte.YEAR] = [
+        sum(hp.energy_consumption[cte.HEATING][cte.MONTH])]
+      boiler.energy_consumption[cte.HEATING][cte.MONTH] = MonthlyValues.get_total_month(
+        boiler.energy_consumption[cte.HEATING][cte.HOUR])
+      boiler.energy_consumption[cte.HEATING][cte.YEAR] = [
+        sum(boiler.energy_consumption[cte.HEATING][cte.MONTH])]
+    if tes.heating_coil_capacity is not None:
+      tes.heating_coil_energy_consumption[cte.HEATING] = {}
+      if len(self.heating_demand) == 8760:
+        tes.heating_coil_energy_consumption[cte.HEATING][cte.HOUR] = coil_hourly
+        tes.heating_coil_energy_consumption[cte.HEATING][cte.MONTH] = MonthlyValues.get_total_month(
+          tes.heating_coil_energy_consumption[cte.HEATING][cte.HOUR])
+        tes.heating_coil_energy_consumption[cte.HEATING][cte.YEAR] = [
+          sum(tes.heating_coil_energy_consumption[cte.HEATING][cte.MONTH])]
+    self.results['Heating Demand (W)'] = demand
+    self.results['HP Heat Output (W)'] = q_hp
+    self.results['HP Source Temperature'] = source_temperature
+    self.results['HP Supply Temperature'] = t_sup_hp
+    self.results['HP COP'] = hp_cop
+    self.results['HP Electricity Consumption (W)'] = hp_electricity
+    self.results['Boiler Heat Output (W)'] = q_boiler
+    self.results['Boiler Power Consumption (W)'] = boiler_energy_consumption
+    self.results['Boiler Supply Temperature'] = t_sup_boiler
+    self.results['Boiler Fuel Consumption'] = boiler_fuel_consumption
+    self.results['TES Temperature'] = t_tank
+    self.results['Heating Coil heat input'] = q_coil
+    self.results['TES Charging Flow Rate (kg/s)'] = m_ch
+    self.results['TES Discharge Flow Rate (kg/s)'] = m_dis
+    self.results['Heating Loop Return Temperature'] = t_ret
+    self.results['Total Heating Power Consumption (W)'] = total_consumption
+    return self.results

energy_system_modelling_package/energy_system_modelling_factories/hvac_dhw_systems_simulation_models/heat_pump_characteristics.py

@@ -0,0 +1,54 @@
+import hub.helpers.constants as cte
+
+
+class HeatPump:
+  def __init__(self, hp, t_out):
+    self.hp = hp
+    self.t_out = t_out
+
+  def hp_source_temperature(self):
+    if self.hp.source_medium == cte.AIR:
+      self.hp.source_temperature = self.t_out
+    elif self.hp.source_medium == cte.GROUND:
+      average_air_temperature = sum(self.t_out) / len(self.t_out)
+      self.hp.source_temperature = [average_air_temperature + 10] * len(self.t_out)
+    elif self.hp.source_medium == cte.WATER:
+      self.hp.source_temperature = [15] * len(self.t_out)
+    return self.hp.source_temperature
+
+  def air_to_water_cop(self, source_temperature, inlet_water_temperature, mode=cte.HEATING):
+    cop_coefficient = 1
+    t_inlet_water_fahrenheit = 1.8 * inlet_water_temperature + 32
+    t_source_fahrenheit = 1.8 * source_temperature + 32
+    if mode == cte.HEATING:
+      if self.hp.heat_efficiency_curve is not None:
+        cop_curve_coefficients = [float(coefficient) for coefficient in self.hp.heat_efficiency_curve.coefficients]
+        cop_coefficient = (1 / (cop_curve_coefficients[0] +
+                                cop_curve_coefficients[1] * t_inlet_water_fahrenheit +
+                                cop_curve_coefficients[2] * t_inlet_water_fahrenheit ** 2 +
+                                cop_curve_coefficients[3] * t_source_fahrenheit +
+                                cop_curve_coefficients[4] * t_source_fahrenheit ** 2 +
+                                cop_curve_coefficients[5] * t_inlet_water_fahrenheit * t_source_fahrenheit))
+      hp_efficiency = float(self.hp.heat_efficiency)
+    elif mode == cte.COOLING:
+      if self.hp.cooling_efficiency_curve is not None:
+        cop_curve_coefficients = [float(coefficient) for coefficient in self.hp.cooling_efficiency_curve.coefficients]
+        cop_coefficient = (1 / (cop_curve_coefficients[0] +
+                                cop_curve_coefficients[1] * t_inlet_water_fahrenheit +
+                                cop_curve_coefficients[2] * t_inlet_water_fahrenheit ** 2 +
+                                cop_curve_coefficients[3] * t_source_fahrenheit +
+                                cop_curve_coefficients[4] * t_source_fahrenheit ** 2 +
+                                cop_curve_coefficients[5] * t_inlet_water_fahrenheit * t_source_fahrenheit)) / 3.41214
+      hp_efficiency = float(self.hp.cooling_efficiency)
+    else:
+      if self.hp.heat_efficiency_curve is not None:
+        cop_curve_coefficients = [float(coefficient) for coefficient in self.hp.heat_efficiency_curve.coefficients]
+        cop_coefficient = (cop_curve_coefficients[0] +
+                           cop_curve_coefficients[1] * source_temperature +
+                           cop_curve_coefficients[2] * source_temperature ** 2 +
+                           cop_curve_coefficients[3] * inlet_water_temperature +
+                           cop_curve_coefficients[4] * inlet_water_temperature ** 2 +
+                           cop_curve_coefficients[5] * inlet_water_temperature * source_temperature)
+      hp_efficiency = float(self.hp.heat_efficiency)
+    hp_cop = cop_coefficient * hp_efficiency
+    return hp_cop

energy_system_modelling_package/energy_system_modelling_factories/hvac_dhw_systems_simulation_models/heat_pump_cooling.py

@@ -0,0 +1,89 @@
+import hub.helpers.constants as cte
+from hub.helpers.monthly_values import MonthlyValues
+from energy_system_modelling_package.energy_system_modelling_factories.hvac_dhw_systems_simulation_models.heat_pump_characteristics import HeatPump
+
+
+class HeatPumpCooling:
+  def __init__(self, hp, hourly_cooling_demand_joules, cooling_peak_load_watts, cutoff_temperature, outdoor_temperature,
+               dt=900):
+    self.hp = hp
+    self.cooling_demand = [demand / cte.WATTS_HOUR_TO_JULES for demand in hourly_cooling_demand_joules]
+    self.cooling_peak_load = cooling_peak_load_watts
+    self.cutoff_temperature = cutoff_temperature
+    self.dt = dt
+    self.results = {}
+    self.heat_pump_characteristics = HeatPump(self.hp, outdoor_temperature)
+
+  def simulation(self):
+    source_temperature_hourly = self.heat_pump_characteristics.hp_source_temperature()
+    cooling_efficiency = float(self.hp.cooling_efficiency)
+    number_of_ts = int(cte.HOUR_TO_SECONDS / self.dt)
+    demand = [0] + [x for x in self.cooling_demand for _ in range(number_of_ts)]
+    source_temperature = [0] + [x for x in source_temperature_hourly for _ in range(number_of_ts)]
+    variable_names = ["t_sup_hp", "t_ret", "m", "q_hp", "hp_electricity", "hp_cop"]
+    num_hours = len(demand)
+    variables = {name: [0] * num_hours for name in variable_names}
+    (t_sup_hp, t_ret, m, q_hp, hp_electricity, hp_cop) = [variables[name] for name in variable_names]
+    t_ret[0] = self.cutoff_temperature
+
+    for i in range(1, len(demand)):
+      if demand[i] > 0.15 * self.cooling_peak_load:
+        m[i] = self.hp.nominal_cooling_output / (cte.WATER_HEAT_CAPACITY * 5)
+        if t_ret[i - 1] >= self.cutoff_temperature:
+          if demand[i] < 0.25 * self.cooling_peak_load:
+            q_hp[i] = 0.25 * self.hp.nominal_cooling_output
+          elif demand[i] < 0.5 * self.cooling_peak_load:
+            q_hp[i] = 0.5 * self.hp.nominal_cooling_output
+          else:
+            q_hp[i] = self.hp.nominal_cooling_output
+          t_sup_hp[i] = t_ret[i - 1] - q_hp[i] / (m[i] * cte.WATER_HEAT_CAPACITY)
+        else:
+          q_hp[i] = 0
+          t_sup_hp[i] = t_ret[i - 1]
+        if m[i] == 0:
+          t_ret[i] = t_sup_hp[i]
+        else:
+          t_ret[i] = t_sup_hp[i] + demand[i] / (m[i] * cte.WATER_HEAT_CAPACITY)
+      else:
+        m[i] = 0
+        q_hp[i] = 0
+        t_sup_hp[i] = t_ret[i - 1]
+        t_ret[i] = t_ret[i - 1]
+      if q_hp[i] > 0:
+        if self.hp.source_medium == cte.AIR and self.hp.supply_medium == cte.WATER:
+          hp_cop[i] = self.heat_pump_characteristics.air_to_water_cop(source_temperature[i], t_ret[i],
+                                                                      mode=cte.COOLING)
+          hp_electricity[i] = q_hp[i] / hp_cop[i]
+      else:
+        hp_cop[i] = 0
+        hp_electricity[i] = 0
+    hp_electricity_j = [(x * cte.WATTS_HOUR_TO_JULES) / number_of_ts for x in hp_electricity]
+    hp_hourly = []
+    hp_supply_temperature_hourly = []
+    hp_sum = 0
+    for i in range(1, len(demand)):
+      hp_sum += hp_electricity_j[i]
+      if (i - 1) % number_of_ts == 0:
+        hp_hourly.append(hp_sum)
+        hp_supply_temperature_hourly.append(t_sup_hp[i])
+        hp_sum = 0
+    self.hp.cooling_supply_temperature = hp_supply_temperature_hourly
+    self.hp.energy_consumption[cte.COOLING] = {}
+    self.hp.energy_consumption[cte.COOLING][cte.HOUR] = hp_hourly
+    self.hp.energy_consumption[cte.COOLING][cte.MONTH] = MonthlyValues.get_total_month(
+      self.hp.energy_consumption[cte.COOLING][cte.HOUR])
+    self.hp.energy_consumption[cte.COOLING][cte.YEAR] = [
+      sum(self.hp.energy_consumption[cte.COOLING][cte.MONTH])]
+    self.results['Cooling Demand (W)'] = demand
+    self.results['HP Cooling Output (W)'] = q_hp
+    self.results['HP Cooling Source Temperature'] = source_temperature
+    self.results['HP Cooling Supply Temperature'] = t_sup_hp
+    self.results['HP Cooling COP'] = hp_cop
+    self.results['HP Electricity Consumption'] = hp_electricity
+    self.results['Cooling Loop Flow Rate (kg/s)'] = m
+    self.results['Cooling Loop Return Temperature'] = t_ret
+    return self.results
+
+
+
+

energy_system_modelling_package/energy_system_modelling_factories/hvac_dhw_systems_simulation_models/thermal_storage_tank.py

@@ -0,0 +1,47 @@
+import math
+
+import hub.helpers.constants as cte
+
+
+class StorageTank:
+  """
+  Calculation of the temperature inside a hot water storage tank in the next time step
+  """
+
+  def __init__(self, volume, height, material_layers, heating_coil_capacity, stratification_layer=1):
+    self.volume = volume
+    self.height = height
+    self.materials = material_layers
+    self.number_of_vertical_layers = stratification_layer
+    self.heating_coil_capacity = heating_coil_capacity
+
+  def heat_loss_coefficient(self):
+    r_tot = sum(float(layer.thickness) / float(layer.material.conductivity) for layer in
+                self.materials)
+    u_tot = 1 / r_tot
+    d = math.sqrt((4 * self.volume) / (math.pi * self.height))
+    a_side = math.pi * d * self.height
+    a_top = math.pi * d ** 2 / 4
+    if self.number_of_vertical_layers == 1:
+      ua = u_tot * (2 * a_top + a_side)
+      return ua
+    else:
+      ua_side = u_tot * a_side
+      ua_top_bottom = u_tot * (a_top + a_side)
+      return ua_side, ua_top_bottom
+
+
+  def calculate_space_heating_fully_mixed(self, charging_flow_rate, discharging_flow_rate, supply_temperature,
+                                          return_temperature,
+                                          current_tank_temperature, heat_generator_input, ambient_temperature, dt):
+    ua = self.heat_loss_coefficient()
+    t_tank = (current_tank_temperature +
+              (charging_flow_rate * (supply_temperature - current_tank_temperature) +
+               (ua * (ambient_temperature - current_tank_temperature)) / cte.WATER_HEAT_CAPACITY -
+               discharging_flow_rate * (current_tank_temperature - return_temperature) +
+               heat_generator_input / cte.WATER_HEAT_CAPACITY) * (dt / (cte.WATER_DENSITY * self.volume)))
+    return t_tank
+
+  def calculate_dhw_fully_mixed(self, charging_flow_rate, discharging_flow_rate, supply_temperature, return_temperature,
+                                current_tank_temperature, heat_generator_input, ambient_temperature, dt):
+    pass

energy_system_modelling_package/energy_system_modelling_factories/montreal_energy_system_archetype_modelling_factory.py

@@ -0,0 +1,36 @@
+"""
+EnergySystemSizingSimulationFactory retrieve the energy system archetype sizing and simulation module
+SPDX - License - Identifier: LGPL - 3.0 - or -later
+Copyright © 2024 Concordia CERC group
+Project Coder Saeed Ranjbar saeed.ranjbar@mail.concordia.ca
+"""
+
+from energy_system_modelling_package.energy_system_modelling_factories.archetypes.montreal.archetype_cluster_1 import ArchetypeCluster1
+
+
+class MontrealEnergySystemArchetypesSimulationFactory:
+  """
+  EnergySystemsFactory class
+  """
+
+  def __init__(self, handler, building, output_path, csv_output=True):
+    self._output_path = output_path
+    self._handler = '_' + handler.lower()
+    self._building = building
+    self._csv_output = csv_output
+
+  def _archetype_cluster_1(self):
+    """
+    Enrich the city by using the sizing and simulation model developed for archetype13 of montreal_future_systems
+    """
+    dt = 900
+    ArchetypeCluster1(self._building, dt, self._output_path, self._csv_output).enrich_building()
+    self._building.level_of_detail.energy_systems = 2
+    self._building.level_of_detail.energy_systems = 2
+
+  def enrich(self):
+    """
+    Enrich the city given to the class using the class given handler
+    :return: None
+    """
+    getattr(self, self._handler, lambda: None)()

energy_system_modelling_package/energy_system_modelling_factories/pv_assessment/electricity_demand_calculator.py

@@ -0,0 +1,73 @@
+import hub.helpers.constants as cte
+class HourlyElectricityDemand:
+  def __init__(self, building):
+    self.building = building
+
+  def calculate(self):
+    hourly_electricity_consumption = []
+    energy_systems = self.building.energy_systems
+    appliance = self.building.appliances_electrical_demand[cte.HOUR]
+    lighting = self.building.lighting_electrical_demand[cte.HOUR]
+    elec_heating = 0
+    elec_cooling = 0
+    elec_dhw = 0
+    if cte.HEATING in self.building.energy_consumption_breakdown[cte.ELECTRICITY]:
+      elec_heating = 1
+    if cte.COOLING in self.building.energy_consumption_breakdown[cte.ELECTRICITY]:
+      elec_cooling = 1
+    if cte.DOMESTIC_HOT_WATER in self.building.energy_consumption_breakdown[cte.ELECTRICITY]:
+      elec_dhw = 1
+    heating = None
+    cooling = None
+    dhw = None
+
+    if elec_heating == 1:
+      for energy_system in energy_systems:
+        if cte.HEATING in energy_system.demand_types:
+          for generation_system in energy_system.generation_systems:
+            if generation_system.fuel_type == cte.ELECTRICITY:
+              if cte.HEATING in generation_system.energy_consumption:
+                heating = generation_system.energy_consumption[cte.HEATING][cte.HOUR]
+              else:
+                if len(energy_system.generation_systems) > 1:
+                  heating = [x / 2 for x in self.building.heating_consumption[cte.HOUR]]
+                else:
+                  heating = self.building.heating_consumption[cte.HOUR]
+
+    if elec_dhw == 1:
+      for energy_system in energy_systems:
+        if cte.DOMESTIC_HOT_WATER in energy_system.demand_types:
+          for generation_system in energy_system.generation_systems:
+            if generation_system.fuel_type == cte.ELECTRICITY:
+              if cte.DOMESTIC_HOT_WATER in generation_system.energy_consumption:
+                dhw = generation_system.energy_consumption[cte.DOMESTIC_HOT_WATER][cte.HOUR]
+              else:
+                if len(energy_system.generation_systems) > 1:
+                  dhw = [x / 2 for x in self.building.domestic_hot_water_consumption[cte.HOUR]]
+                else:
+                  dhw = self.building.domestic_hot_water_consumption[cte.HOUR]
+
+    if elec_cooling == 1:
+      for energy_system in energy_systems:
+        if cte.COOLING in energy_system.demand_types:
+          for generation_system in energy_system.generation_systems:
+            if cte.COOLING in generation_system.energy_consumption:
+              cooling = generation_system.energy_consumption[cte.COOLING][cte.HOUR]
+            else:
+              if len(energy_system.generation_systems) > 1:
+                cooling = [x / 2 for x in self.building.cooling_consumption[cte.HOUR]]
+              else:
+                cooling = self.building.cooling_consumption[cte.HOUR]
+
+    for i in range(len(self.building.heating_demand[cte.HOUR])):
+      hourly = 0
+      hourly += appliance[i]
+      hourly += lighting[i]
+      if heating is not None:
+        hourly += heating[i]
+      if cooling is not None:
+        hourly += cooling[i]
+      if dhw is not None:
+        hourly += dhw[i]
+      hourly_electricity_consumption.append(hourly)
+    return hourly_electricity_consumption

energy_system_modelling_package/energy_system_modelling_factories/pv_assessment/pv_feasibility.py

@@ -0,0 +1,37 @@
+from pathlib import Path
+import subprocess
+from hub.imports.geometry_factory import GeometryFactory
+from building_modelling.geojson_creator import process_geojson
+from hub.helpers.dictionaries import Dictionaries
+from hub.imports.weather_factory import WeatherFactory
+from hub.imports.results_factory import ResultFactory
+from hub.exports.exports_factory import ExportsFactory
+
+
+def pv_feasibility(current_x, current_y, current_diff, selected_buildings):
+  input_files_path = (Path(__file__).parent.parent.parent.parent / 'input_files')
+  output_path = (Path(__file__).parent.parent.parent.parent / 'out_files').resolve()
+  sra_output_path = output_path / 'sra_outputs' / 'extended_city_sra_outputs'
+  sra_output_path.mkdir(parents=True, exist_ok=True)
+  new_diff = current_diff * 5
+  geojson_file = process_geojson(x=current_x, y=current_y, diff=new_diff, expansion=True)
+  file_path = input_files_path / 'output_buildings.geojson'
+  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
+  WeatherFactory('epw', city).enrich()
+  ExportsFactory('sra', city, sra_output_path).export()
+  sra_path = (sra_output_path / f'{city.name}_sra.xml').resolve()
+  subprocess.run(['sra', str(sra_path)])
+  ResultFactory('sra', city, sra_output_path).enrich()
+  for selected_building in selected_buildings:
+    for building in city.buildings:
+      if selected_building.name == building.name:
+        selected_building.roofs[0].global_irradiance = building.roofs[0].global_irradiance
+
+
+
+

energy_system_modelling_package/energy_system_modelling_factories/pv_assessment/pv_model.py

@@ -0,0 +1,42 @@
+import math
+import hub.helpers.constants as cte
+from hub.helpers.monthly_values import MonthlyValues
+
+
+class PVModel:
+  def __init__(self, pv, hourly_electricity_demand_joules, solar_radiation, installed_pv_area, model_type, ns=None,
+               np=None):
+    self.pv = pv
+    self.hourly_electricity_demand = [demand / cte.WATTS_HOUR_TO_JULES for demand in hourly_electricity_demand_joules]
+    self.solar_radiation = solar_radiation
+    self.installed_pv_area = installed_pv_area
+    self._model_type = '_' + model_type.lower()
+    self.ns = ns
+    self.np = np
+    self.results = {}
+
+  def fixed_efficiency(self):
+    module_efficiency = float(self.pv.electricity_efficiency)
+    variable_names = ["pv_output", "import", "export", "self_sufficiency_ratio"]
+    variables = {name: [0] * len(self.hourly_electricity_demand) for name in variable_names}
+    (pv_out, grid_import, grid_export, self_sufficiency_ratio) = [variables[name] for name in variable_names]
+    for i in range(len(self.hourly_electricity_demand)):
+      pv_out[i] = module_efficiency * self.installed_pv_area * self.solar_radiation[i] / cte.WATTS_HOUR_TO_JULES
+      if pv_out[i] < self.hourly_electricity_demand[i]:
+        grid_import[i] = self.hourly_electricity_demand[i] - pv_out[i]
+      else:
+        grid_export[i] = pv_out[i] - self.hourly_electricity_demand[i]
+      self_sufficiency_ratio[i] = pv_out[i] / self.hourly_electricity_demand[i]
+    self.results['Electricity Demand (W)'] = self.hourly_electricity_demand
+    self.results['PV Output (W)'] = pv_out
+    self.results['Imported from Grid (W)'] = grid_import
+    self.results['Exported to Grid (W)'] = grid_export
+    self.results['Self Sufficiency Ratio'] = self_sufficiency_ratio
+    return self.results
+
+  def enrich(self):
+    """
+    Enrich the city given to the class using the class given handler
+    :return: None
+    """
+    return getattr(self, self._model_type, lambda: None)()

energy_system_modelling_package/energy_system_modelling_factories/pv_assessment/pv_sizing.py

@@ -0,0 +1,70 @@
+import math
+import hub.helpers.constants as cte
+from energy_system_modelling_package.energy_system_modelling_factories.pv_assessment.solar_angles import CitySolarAngles
+from energy_system_modelling_package.energy_system_modelling_factories.pv_assessment.radiation_tilted import RadiationTilted
+
+
+class PVSizing(CitySolarAngles):
+  def __init__(self, city, tilt_angle, surface_azimuth=180, maintenance_factor=0.1, mechanical_equipment_factor=0.3,
+               orientation_factor=0.1, system_type='rooftop'):
+    super().__init__(location_latitude=city.latitude,
+                     location_longitude=city.longitude,
+                     tilt_angle=tilt_angle,
+                     surface_azimuth_angle=surface_azimuth)
+    self.city = city
+    self.maintenance_factor = maintenance_factor
+    self.mechanical_equipment_factor = mechanical_equipment_factor
+    self.orientation_factor = orientation_factor
+    self.angles = self.calculate
+    self.system_type = system_type
+
+  def rooftop_sizing(self):
+    results = {}
+    # Available Roof Area
+    for building in self.city.buildings:
+      for energy_system in building.energy_systems:
+        for generation_system in energy_system.generation_systems:
+          if generation_system.system_type == cte.PHOTOVOLTAIC:
+            module_width = float(generation_system.width)
+            module_height = float(generation_system.height)
+            roof_area = 0
+            for roof in building.roofs:
+              roof_area += roof.perimeter_area
+            pv_module_area = module_width * module_height
+            available_roof = ((self.maintenance_factor + self.orientation_factor + self.mechanical_equipment_factor) *
+                              roof_area)
+        # Inter-Row Spacing
+            winter_solstice = self.angles[(self.angles['AST'].dt.month == 12) &
+                                          (self.angles['AST'].dt.day == 21) &
+                                          (self.angles['AST'].dt.hour == 12)]
+            solar_altitude = winter_solstice['solar altitude'].values[0]
+            solar_azimuth = winter_solstice['solar azimuth'].values[0]
+            distance = ((module_height * abs(math.cos(math.radians(solar_azimuth)))) /
+                        math.tan(math.radians(solar_altitude)))
+            distance = float(format(distance, '.1f'))
+        # Calculation of the number of panels
+            space_dimension = math.sqrt(available_roof)
+            space_dimension = float(format(space_dimension, '.2f'))
+            panels_per_row = math.ceil(space_dimension / module_width)
+            number_of_rows = math.ceil(space_dimension / distance)
+            total_number_of_panels = panels_per_row * number_of_rows
+            total_pv_area = panels_per_row * number_of_rows * pv_module_area
+            building.roofs[0].installed_solar_collector_area = total_pv_area
+            results[f'Building {building.name}'] = {'total_roof_area': roof_area,
+                                                    'PV dedicated area': available_roof,
+                                                    'total_pv_area': total_pv_area,
+                                                    'total_number_of_panels': total_number_of_panels,
+                                                    'number_of_rows': number_of_rows,
+                                                    'panels_per_row': panels_per_row}
+    return results
+
+  def rooftop_tilted_radiation(self):
+    for building in self.city.buildings:
+      RadiationTilted(building=building,
+                      solar_angles=self.angles,
+                      tilt_angle=self.tilt_angle,
+                      ghi=building.roofs[0].global_irradiance[cte.HOUR],
+                      ).enrich()
+
+  def facade_sizing(self):
+    pass

energy_system_modelling_package/energy_system_modelling_factories/pv_assessment/pv_sizing_and_simulation.py

@@ -0,0 +1,59 @@
+import math
+
+from energy_system_modelling_package.energy_system_modelling_factories.pv_assessment.radiation_tilted import RadiationTilted
+import hub.helpers.constants as cte
+from hub.helpers.monthly_values import MonthlyValues
+
+
+class PVSizingSimulation(RadiationTilted):
+  def __init__(self, building, solar_angles, tilt_angle, module_height, module_width, ghi):
+    super().__init__(building, solar_angles, tilt_angle, ghi)
+    self.module_height = module_height
+    self.module_width = module_width
+    self.total_number_of_panels = 0
+    self.enrich()
+
+  def available_space(self):
+    roof_area = self.building.roofs[0].perimeter_area
+    maintenance_factor = 0.1
+    orientation_factor = 0.2
+    if self.building.function == cte.RESIDENTIAL:
+      mechanical_equipment_factor = 0.2
+    else:
+      mechanical_equipment_factor = 0.3
+    available_roof = (maintenance_factor + orientation_factor + mechanical_equipment_factor) * roof_area
+    return available_roof
+
+  def inter_row_spacing(self):
+    winter_solstice = self.df[(self.df['AST'].dt.month == 12) &
+                              (self.df['AST'].dt.day == 21) &
+                              (self.df['AST'].dt.hour == 12)]
+    solar_altitude = winter_solstice['solar altitude'].values[0]
+    solar_azimuth = winter_solstice['solar azimuth'].values[0]
+    distance = ((self.module_height * abs(math.cos(math.radians(solar_azimuth)))) /
+                math.tan(math.radians(solar_altitude)))
+    distance = float(format(distance, '.1f'))
+    return distance
+
+  def number_of_panels(self, available_roof, inter_row_distance):
+    space_dimension = math.sqrt(available_roof)
+    space_dimension = float(format(space_dimension, '.2f'))
+    panels_per_row = math.ceil(space_dimension / self.module_width)
+    number_of_rows = math.ceil(space_dimension / inter_row_distance)
+    self.total_number_of_panels = panels_per_row * number_of_rows
+    return panels_per_row, number_of_rows
+
+  def pv_output_constant_efficiency(self):
+    radiation = self.total_radiation_tilted
+    pv_module_area = self.module_width * self.module_height
+    available_roof = self.available_space()
+    inter_row_spacing = self.inter_row_spacing()
+    self.number_of_panels(available_roof, inter_row_spacing)
+    self.building.roofs[0].installed_solar_collector_area = pv_module_area * self.total_number_of_panels
+    system_efficiency = 0.2
+    pv_hourly_production = [x * system_efficiency * self.total_number_of_panels * pv_module_area *
+                            cte.WATTS_HOUR_TO_JULES for x in radiation]
+    self.building.onsite_electrical_production[cte.HOUR] = pv_hourly_production
+    self.building.onsite_electrical_production[cte.MONTH] = (
+      MonthlyValues.get_total_month(self.building.onsite_electrical_production[cte.HOUR]))
+    self.building.onsite_electrical_production[cte.YEAR] = [sum(self.building.onsite_electrical_production[cte.MONTH])]

energy_system_modelling_package/energy_system_modelling_factories/pv_assessment/radiation_tilted.py

@@ -0,0 +1,110 @@
+import pandas as pd
+import math
+import hub.helpers.constants as cte
+from hub.helpers.monthly_values import MonthlyValues
+
+
+class RadiationTilted:
+  def __init__(self, building, solar_angles, tilt_angle, ghi, solar_constant=1366.1, maximum_clearness_index=1,
+               min_cos_zenith=0.065, maximum_zenith_angle=87):
+    self.building = building
+    self.ghi = ghi
+    self.tilt_angle = tilt_angle
+    self.zeniths = solar_angles['zenith'].tolist()[:-1]
+    self.incidents = solar_angles['incident angle'].tolist()[:-1]
+    self.date_time = solar_angles['DateTime'].tolist()[:-1]
+    self.ast = solar_angles['AST'].tolist()[:-1]
+    self.solar_azimuth = solar_angles['solar azimuth'].tolist()[:-1]
+    self.solar_altitude = solar_angles['solar altitude'].tolist()[:-1]
+    data = {'DateTime': self.date_time, 'AST': self.ast, 'solar altitude': self.solar_altitude, 'zenith': self.zeniths,
+            'solar azimuth': self.solar_azimuth, 'incident angle': self.incidents, 'ghi': self.ghi}
+    self.df = pd.DataFrame(data)
+    self.df['DateTime'] = pd.to_datetime(self.df['DateTime'])
+    self.df['AST'] = pd.to_datetime(self.df['AST'])
+    self.df.set_index('DateTime', inplace=True)
+    self.solar_constant = solar_constant
+    self.maximum_clearness_index = maximum_clearness_index
+    self.min_cos_zenith = min_cos_zenith
+    self.maximum_zenith_angle = maximum_zenith_angle
+    self.i_on = []
+    self.i_oh = []
+    self.k_t = []
+    self.fraction_diffuse = []
+    self.diffuse_horizontal = []
+    self.beam_horizontal = []
+    self.dni = []
+    self.beam_tilted = []
+    self.diffuse_tilted = []
+    self.total_radiation_tilted = []
+    self.calculate()
+
+  def dni_extra(self):
+    for i in range(len(self.df)):
+      self.i_on.append(self.solar_constant * (1 + 0.033 * math.cos(math.radians(360 * self.df.index.dayofyear[i] / 365))))
+
+    self.df['extraterrestrial normal radiation (Wh/m2)'] = self.i_on
+
+  def clearness_index(self):
+    for i in range(len(self.df)):
+      self.i_oh.append(self.i_on[i] * max(math.cos(math.radians(self.zeniths[i])), self.min_cos_zenith))
+      self.k_t.append(self.ghi[i] / self.i_oh[i])
+      self.k_t[i] = max(0, self.k_t[i])
+      self.k_t[i] = min(self.maximum_clearness_index, self.k_t[i])
+    self.df['extraterrestrial radiation on horizontal (Wh/m2)'] = self.i_oh
+    self.df['clearness index'] = self.k_t
+
+  def diffuse_fraction(self):
+    for i in range(len(self.df)):
+      if self.k_t[i] <= 0.22:
+        self.fraction_diffuse.append(1 - 0.09 * self.k_t[i])
+      elif self.k_t[i] <= 0.8:
+        self.fraction_diffuse.append(0.9511 - 0.1604 * self.k_t[i] + 4.388 * self.k_t[i] ** 2 -
+                                     16.638 * self.k_t[i] ** 3 + 12.336 * self.k_t[i] ** 4)
+      else:
+        self.fraction_diffuse.append(0.165)
+      if self.zeniths[i] > self.maximum_zenith_angle:
+        self.fraction_diffuse[i] = 1
+
+    self.df['diffuse fraction'] = self.fraction_diffuse
+
+  def radiation_components_horizontal(self):
+    for i in range(len(self.df)):
+      self.diffuse_horizontal.append(self.ghi[i] * self.fraction_diffuse[i])
+      self.beam_horizontal.append(self.ghi[i] - self.diffuse_horizontal[i])
+      self.dni.append((self.ghi[i] - self.diffuse_horizontal[i]) / math.cos(math.radians(self.zeniths[i])))
+      if self.zeniths[i] > self.maximum_zenith_angle or self.dni[i] < 0:
+        self.dni[i] = 0
+
+    self.df['diffuse horizontal (Wh/m2)'] = self.diffuse_horizontal
+    self.df['dni (Wh/m2)'] = self.dni
+    self.df['beam horizontal (Wh/m2)'] = self.beam_horizontal
+
+  def radiation_components_tilted(self):
+    for i in range(len(self.df)):
+      self.beam_tilted.append(self.dni[i] * math.cos(math.radians(self.incidents[i])))
+      self.beam_tilted[i] = max(self.beam_tilted[i], 0)
+      self.diffuse_tilted.append(self.diffuse_horizontal[i] * ((1 + math.cos(math.radians(self.tilt_angle))) / 2))
+      self.total_radiation_tilted.append(self.beam_tilted[i] + self.diffuse_tilted[i])
+
+    self.df['beam tilted (Wh/m2)'] = self.beam_tilted
+    self.df['diffuse tilted (Wh/m2)'] = self.diffuse_tilted
+    self.df['total radiation tilted (Wh/m2)'] = self.total_radiation_tilted
+
+  def calculate(self) -> pd.DataFrame:
+    self.dni_extra()
+    self.clearness_index()
+    self.diffuse_fraction()
+    self.radiation_components_horizontal()
+    self.radiation_components_tilted()
+    return self.df
+
+  def enrich(self):
+    tilted_radiation = self.total_radiation_tilted
+    self.building.roofs[0].global_irradiance_tilted[cte.HOUR] = tilted_radiation
+    self.building.roofs[0].global_irradiance_tilted[cte.MONTH] = (
+      MonthlyValues.get_total_month(self.building.roofs[0].global_irradiance_tilted[cte.HOUR]))
+    self.building.roofs[0].global_irradiance_tilted[cte.YEAR] = \
+        [sum(self.building.roofs[0].global_irradiance_tilted[cte.MONTH])]
+
+
+

energy_system_modelling_package/energy_system_modelling_factories/pv_assessment/solar_angles.py

@@ -0,0 +1,145 @@
+import math
+import pandas as pd
+from datetime import datetime
+from pathlib import Path
+
+
+class CitySolarAngles:
+  def __init__(self, location_latitude, location_longitude, tilt_angle, surface_azimuth_angle,
+               standard_meridian=-75):
+    self.location_latitude = location_latitude
+    self.location_longitude = location_longitude
+    self.location_latitude_rad = math.radians(location_latitude)
+    self.surface_azimuth_angle = surface_azimuth_angle
+    self.surface_azimuth_rad = math.radians(surface_azimuth_angle)
+    self.tilt_angle = tilt_angle
+    self.tilt_angle_rad = math.radians(tilt_angle)
+    self.standard_meridian = standard_meridian
+    self.longitude_correction = (location_longitude - standard_meridian) * 4
+    self.timezone = 'Etc/GMT+5'
+
+    self.eot = []
+    self.ast = []
+    self.hour_angles = []
+    self.declinations = []
+    self.solar_altitudes = []
+    self.solar_azimuths = []
+    self.zeniths = []
+    self.incidents = []
+    self.beam_tilted = []
+    self.factor = []
+    self.times = pd.date_range(start='2023-01-01', end='2024-01-01', freq='h', tz=self.timezone)
+    self.df = pd.DataFrame(index=self.times)
+    self.day_of_year = self.df.index.dayofyear
+
+  def solar_time(self, datetime_val, day_of_year):
+    b = (day_of_year - 81) * 2 * math.pi / 364
+    eot = 9.87 * math.sin(2 * b) - 7.53 * math.cos(b) - 1.5 * math.sin(b)
+    self.eot.append(eot)
+
+    # Calculate Local Solar Time (LST)
+    lst_hour = datetime_val.hour
+    lst_minute = datetime_val.minute
+    lst_second = datetime_val.second
+    lst = lst_hour + lst_minute / 60 + lst_second / 3600
+
+    # Calculate Apparent Solar Time (AST) in decimal hours
+    ast_decimal = lst + eot / 60 + self.longitude_correction / 60
+    ast_hours = int(ast_decimal)
+    ast_minutes = round((ast_decimal - ast_hours) * 60)
+
+    # Ensure ast_minutes is within valid range
+    if ast_minutes == 60:
+      ast_hours += 1
+      ast_minutes = 0
+    elif ast_minutes < 0:
+      ast_minutes = 0
+    ast_time = datetime(year=datetime_val.year, month=datetime_val.month, day=datetime_val.day,
+                        hour=ast_hours, minute=ast_minutes)
+    self.ast.append(ast_time)
+    return ast_time
+
+  def declination_angle(self, day_of_year):
+    declination = 23.45 * math.sin(math.radians(360 / 365 * (284 + day_of_year)))
+    declination_radian = math.radians(declination)
+    self.declinations.append(declination)
+    return declination_radian
+
+  def hour_angle(self, ast_time):
+    hour_angle = ((ast_time.hour * 60 + ast_time.minute) - 720) / 4
+    hour_angle_radian = math.radians(hour_angle)
+    self.hour_angles.append(hour_angle)
+    return hour_angle_radian
+
+  def solar_altitude(self, declination_radian, hour_angle_radian):
+    solar_altitude_radians = math.asin(math.cos(self.location_latitude_rad) * math.cos(declination_radian) *
+                                       math.cos(hour_angle_radian) + math.sin(self.location_latitude_rad) *
+                                       math.sin(declination_radian))
+    solar_altitude = math.degrees(solar_altitude_radians)
+    self.solar_altitudes.append(solar_altitude)
+    return solar_altitude_radians
+
+  def zenith(self, solar_altitude_radians):
+    solar_altitude = math.degrees(solar_altitude_radians)
+    zenith_degree = 90 - solar_altitude
+    zenith_radian = math.radians(zenith_degree)
+    self.zeniths.append(zenith_degree)
+    return zenith_radian
+
+  def solar_azimuth_analytical(self, hourangle, declination, zenith):
+    numer = (math.cos(zenith) * math.sin(self.location_latitude_rad) - math.sin(declination))
+    denom = (math.sin(zenith) * math.cos(self.location_latitude_rad))
+    if math.isclose(denom, 0.0, abs_tol=1e-8):
+      cos_azi = 1.0
+    else:
+      cos_azi = numer / denom
+
+    cos_azi = max(-1.0, min(1.0, cos_azi))
+
+    sign_ha = math.copysign(1, hourangle)
+    solar_azimuth_radians = sign_ha * math.acos(cos_azi) + math.pi
+    solar_azimuth_degrees = math.degrees(solar_azimuth_radians)
+    self.solar_azimuths.append(solar_azimuth_degrees)
+    return solar_azimuth_radians
+
+  def incident_angle(self, solar_altitude_radians, solar_azimuth_radians):
+    incident_radian = math.acos(math.cos(solar_altitude_radians) *
+                                math.cos(abs(solar_azimuth_radians - self.surface_azimuth_rad)) *
+                                math.sin(self.tilt_angle_rad) + math.sin(solar_altitude_radians) *
+                                math.cos(self.tilt_angle_rad))
+    incident_angle_degrees = math.degrees(incident_radian)
+    self.incidents.append(incident_angle_degrees)
+    return incident_radian
+
+  @property
+  def calculate(self) -> pd.DataFrame:
+    for i in range(len(self.times)):
+      datetime_val = self.times[i]
+      day_of_year = self.day_of_year[i]
+      declination_radians = self.declination_angle(day_of_year)
+      ast_time = self.solar_time(datetime_val, day_of_year)
+      hour_angle_radians = self.hour_angle(ast_time)
+      solar_altitude_radians = self.solar_altitude(declination_radians, hour_angle_radians)
+      zenith_radians = self.zenith(solar_altitude_radians)
+      solar_azimuth_radians = self.solar_azimuth_analytical(hour_angle_radians, declination_radians, zenith_radians)
+      incident_angle_radian = self.incident_angle(solar_altitude_radians, solar_azimuth_radians)
+
+    self.df['DateTime'] = self.times
+    self.df['AST'] = self.ast
+    self.df['hour angle'] = self.hour_angles
+    self.df['eot'] = self.eot
+    self.df['declination angle'] = self.declinations
+    self.df['solar altitude'] = self.solar_altitudes
+    self.df['zenith'] = self.zeniths
+    self.df['solar azimuth'] = self.solar_azimuths
+    self.df['incident angle'] = self.incidents
+
+    return self.df
+
+
+
+
+
+
+
+

energy_system_modelling_package/energy_system_modelling_factories/system_sizing_methods/genetic_algorithm/individual.py

@@ -0,0 +1,470 @@
+import math
+import random
+from hub.helpers.dictionaries import Dictionaries
+from hub.catalog_factories.costs_catalog_factory import CostsCatalogFactory
+import hub.helpers.constants as cte
+from energy_system_modelling_package.energy_system_modelling_factories.hvac_dhw_systems_simulation_models.domestic_hot_water_heat_pump_with_tes import \
+  DomesticHotWaterHeatPumpTes
+from energy_system_modelling_package.energy_system_modelling_factories.hvac_dhw_systems_simulation_models.heat_pump_boiler_tes_heating import \
+  HeatPumpBoilerTesHeating
+import numpy_financial as npf
+
+
+class Individual:
+  def __init__(self, building, energy_system, design_period_energy_demands, optimization_scenario,
+               heating_design_load=None, cooling_design_load=None, dt=900, fuel_price_index=0.05,
+               electricity_tariff_type='fixed', consumer_price_index=0.04, interest_rate=0.04,
+               discount_rate=0.03, percentage_credit=0, credit_years=15):
+    """
+    :param building: building object
+    :param energy_system: energy system to be optimized
+    :param design_period_energy_demands: A dictionary of design period heating, cooling and dhw demands. Design period
+    is the day with the highest total demand and the two days before and after it
+    :param optimization_scenario: a string indicating the objective function from minimization of cost,
+    energy consumption, and both together
+    :param heating_design_load: heating design load in W
+    :param cooling_design_load: cooling design load in W
+    :param dt the time step size used for simulations
+    :param fuel_price_index the price increase index of all fuels. A single value is used for all fuels.
+    :param electricity_tariff_type the electricity tariff type between 'fixed' and 'variable' for economic optimization
+    :param consumer_price_index
+    """
+    self.building = building
+    self.energy_system = energy_system
+    self.design_period_energy_demands = design_period_energy_demands
+    self.demand_types = energy_system.demand_types
+    self.optimization_scenario = optimization_scenario
+    self.heating_design_load = heating_design_load
+    self.cooling_design_load = cooling_design_load
+    self.available_space = building.volume / building.storeys_above_ground
+    self.dt = dt
+    self.fuel_price_index = fuel_price_index
+    self.electricity_tariff_type = electricity_tariff_type
+    self.consumer_price_index = consumer_price_index
+    self.interest_rate = interest_rate
+    self.discount_rate = discount_rate
+    self.credit_years = credit_years
+    self.percentage_credit = percentage_credit
+    self.costs = self.costs_archetype()
+    self.feasibility = True
+    self.fitness_score = 0
+    self.individual = {}
+
+  def system_components(self):
+    """
+    Extracts system components (generation and storage) for a given energy system.
+    :return: Dictionary of system components.
+    """
+    self.individual['Generation Components'] = []
+    self.individual['Energy Storage Components'] = []
+    self.individual['End of Life Cost'] = self.costs.end_of_life_cost
+    for generation_component in self.energy_system.generation_systems:
+      investment_cost, reposition_cost, lifetime = self.unit_investment_cost('Generation',
+                                                                             generation_component.system_type)
+      maintenance_cost = self.unit_maintenance_cost(generation_component)
+      if generation_component.system_type == cte.PHOTOVOLTAIC:
+        heating_capacity = None
+        cooling_capacity = None
+        heat_efficiency = None
+        cooling_efficiency = None
+        unit_fuel_cost = 0
+      else:
+        heating_capacity = 0
+        cooling_capacity = 0
+        heat_efficiency = generation_component.heat_efficiency
+        cooling_efficiency = generation_component.cooling_efficiency
+        unit_fuel_cost = self.fuel_cost_per_kwh(generation_component.fuel_type, 'fixed')
+      self.individual['Generation Components'].append({
+        'type': generation_component.system_type,
+        'heating_capacity': heating_capacity,
+        'cooling_capacity': cooling_capacity,
+        'electricity_capacity': 0,
+        'nominal_heating_efficiency': heat_efficiency,
+        'nominal_cooling_efficiency': cooling_efficiency,
+        'nominal_electricity_efficiency': generation_component.electricity_efficiency,
+        'fuel_type': generation_component.fuel_type,
+        'unit_investment_cost': investment_cost,
+        'unit_reposition_cost': reposition_cost,
+        'unit_fuel_cost(CAD/kWh)': unit_fuel_cost,
+        'unit_maintenance_cost': maintenance_cost,
+        'lifetime': lifetime
+      })
+      if generation_component.energy_storage_systems is not None:
+        for energy_storage_system in generation_component.energy_storage_systems:
+          investment_cost, reposition_cost, lifetime = (
+            self.unit_investment_cost('Storage',
+                                      f'{energy_storage_system.type_energy_stored}_storage'))
+          if energy_storage_system.type_energy_stored == cte.THERMAL:
+            heating_coil_capacity = energy_storage_system.heating_coil_capacity
+            heating_coil_fuel_cost = self.fuel_cost_per_kwh(f'{cte.ELECTRICITY}', 'fixed')
+            volume = 0
+            capacity = None
+          else:
+            heating_coil_capacity = None
+            heating_coil_fuel_cost = None
+            volume = None
+            capacity = 0
+          self.individual['Energy Storage Components'].append(
+            {'type': f'{energy_storage_system.type_energy_stored}_storage',
+             'capacity': capacity,
+             'volume': volume,
+             'heating_coil_capacity': heating_coil_capacity,
+             'unit_investment_cost': investment_cost,
+             'unit_reposition_cost': reposition_cost,
+             'heating_coil_fuel_cost': heating_coil_fuel_cost,
+             'unit_maintenance_cost': 0,
+             'lifetime': lifetime
+             })
+
+  def initialization(self):
+    """
+    Assigns initial sizes to generation and storage components based on heating and cooling design loads and
+    available space in the building.
+    :return:
+    """
+    self.system_components()
+    generation_components = self.individual['Generation Components']
+    storage_components = self.individual['Energy Storage Components']
+    for generation_component in generation_components:
+      if generation_component[
+        'nominal_heating_efficiency'] is not None and cte.HEATING or cte.DOMESTIC_HOT_WATER in self.demand_types:
+        if self.heating_design_load is not None:
+          generation_component['heating_capacity'] = random.uniform(0, self.heating_design_load)
+        else:
+          if cte.HEATING in self.demand_types:
+            generation_component['heating_capacity'] = random.uniform(0, max(
+              self.design_period_energy_demands[cte.HEATING]['demands']) / cte.WATTS_HOUR_TO_JULES)
+          else:
+            generation_component['heating_capacity'] = random.uniform(0, max(
+              self.design_period_energy_demands[cte.DOMESTIC_HOT_WATER]['demands']) / cte.WATTS_HOUR_TO_JULES)
+      else:
+        generation_component['heating_capacity'] = None
+      if generation_component['nominal_cooling_efficiency'] is not None and cte.COOLING in self.demand_types:
+        if self.cooling_design_load is not None:
+          generation_component['cooling_capacity'] = random.uniform(0, self.cooling_design_load)
+        else:
+          generation_component['cooling_capacity'] = random.uniform(0, max(
+            self.design_period_energy_demands[cte.COOLING]['demands']) / cte.WATTS_HOUR_TO_JULES)
+      else:
+        generation_component['cooling_capacity'] = None
+      if generation_component['nominal_electricity_efficiency'] is None:
+        generation_component['electricity_capacity'] = None
+    for storage_component in storage_components:
+      if storage_component['type'] == f'{cte.THERMAL}_storage':
+        storage_component['volume'] = random.uniform(0, 0.01 * self.available_space)
+        if storage_component['heating_coil_capacity'] is not None:
+          if self.heating_design_load is not None:
+            storage_component['heating_coil_capacity'] = random.uniform(0, self.heating_design_load)
+          else:
+            if cte.HEATING in self.demand_types:
+              storage_component['heating_coil_capacity'] = random.uniform(0, max(
+                self.design_period_energy_demands[cte.HEATING]['demands']) / cte.WATTS_HOUR_TO_JULES)
+            else:
+              storage_component['heating_coil_capacity'] = random.uniform(0, max(
+                self.design_period_energy_demands[cte.DOMESTIC_HOT_WATER]['demands']) / cte.WATTS_HOUR_TO_JULES)
+
+  def score_evaluation(self):
+    self.system_simulation()
+    self.individual['feasible'] = self.feasibility
+    lcc = 0
+    total_energy_consumption = 0
+    if self.feasibility:
+      if 'cost' in self.optimization_scenario:
+        investment_cost = 0
+        operation_cost_year_0 = 0
+        maintenance_cost_year_0 = 0
+        for generation_system in self.individual['Generation Components']:
+          heating_capacity = 0 if generation_system['heating_capacity'] is None else generation_system[
+            'heating_capacity']
+          cooling_capacity = 0 if generation_system['cooling_capacity'] is None else generation_system[
+            'cooling_capacity']
+          capacity = max(heating_capacity, cooling_capacity)
+          investment_cost += capacity * generation_system['unit_investment_cost'] / 1000
+          maintenance_cost_year_0 += capacity * generation_system['unit_maintenance_cost'] / 1000
+          operation_cost_year_0 += (generation_system['total_energy_consumption(kWh)'] *
+                                    generation_system['unit_fuel_cost(CAD/kWh)'])
+        for storage_system in self.individual['Energy Storage Components']:
+          if cte.THERMAL in storage_system['type']:
+            investment_cost += storage_system['volume'] * storage_system['unit_investment_cost']
+            if storage_system['heating_coil_capacity'] is not None:
+              operation_cost_year_0 += (storage_system['total_energy_consumption(kWh)'] *
+                                        storage_system['heating_coil_fuel_cost'])
+        lcc = self.life_cycle_cost_calculation(investment_cost=investment_cost,
+                                               operation_cost_year_0=operation_cost_year_0,
+                                               maintenance_cost_year_0=maintenance_cost_year_0)
+        self.individual['lcc'] = lcc
+      if 'energy_consumption' in self.optimization_scenario:
+        total_energy_consumption = 0
+        for generation_system in self.individual['Generation Components']:
+          total_energy_consumption += generation_system['total_energy_consumption(kWh)']
+        for storage_system in self.individual['Energy Storage Components']:
+          if cte.THERMAL in storage_system['type'] and storage_system['heating_coil_capacity'] is not None:
+            total_energy_consumption += storage_system['total_energy_consumption(kWh)']
+        self.individual['total_energy_consumption'] = total_energy_consumption
+    # Fitness score based on the optimization scenario
+      if self.optimization_scenario == 'cost':
+        self.fitness_score = lcc
+        self.individual['fitness_score'] = lcc
+      elif self.optimization_scenario == 'energy_consumption':
+        self.fitness_score = total_energy_consumption
+        self.individual['fitness_score'] = total_energy_consumption
+      elif self.optimization_scenario == 'cost_energy_consumption' or self.optimization_scenario == 'energy_consumption_cost':
+        self.fitness_score = (lcc, total_energy_consumption)
+        self.individual['fitness_score'] = (lcc, total_energy_consumption)
+    else:
+      lcc = float('inf')
+      total_energy_consumption = float('inf')
+      self.individual['lcc'] = lcc
+      self.individual['total_energy_consumption'] = total_energy_consumption
+      if self.optimization_scenario == 'cost_energy_consumption' or self.optimization_scenario == 'energy_consumption_cost':
+        self.individual['fitness_score'] = (float('inf'), float('inf'))
+        self.fitness_score = (float('inf'), float('inf'))
+      else:
+        self.individual['fitness_score'] = float('inf')
+        self.fitness_score = float('inf')
+
+  def system_simulation(self):
+    """
+    The method to run the energy system model using the existing models in the energy_system_modelling_package.
+    Based on cluster id and demands, model is imported and run.
+    :return: dictionary of results
+    """
+    if self.building.energy_systems_archetype_cluster_id == '1':
+      if cte.HEATING in self.demand_types:
+        boiler = self.energy_system.generation_systems[0]
+        boiler.nominal_heat_output = self.individual['Generation Components'][0]['heating_capacity']
+        hp = self.energy_system.generation_systems[1]
+        hp.nominal_heat_output = self.individual['Generation Components'][1]['heating_capacity']
+        tes = self.energy_system.generation_systems[0].energy_storage_systems[0]
+        tes.volume = self.individual['Energy Storage Components'][0]['volume']
+        tes.height = self.building.average_storey_height - 1
+        tes.heating_coil_capacity = self.individual['Energy Storage Components'][0]['heating_coil_capacity'] \
+          if self.individual['Energy Storage Components'][0]['heating_coil_capacity'] is not None else None
+        heating_demand_joules = self.design_period_energy_demands[cte.HEATING]['demands']
+        heating_peak_load_watts = max(self.design_period_energy_demands[cte.HEATING]) if \
+          (self.heating_design_load is not None) else self.building.heating_peak_load[cte.YEAR][0]
+        upper_limit_tes_heating = 55
+        design_period_start_index = self.design_period_energy_demands[cte.HEATING]['start_index']
+        design_period_end_index = self.design_period_energy_demands[cte.HEATING]['end_index']
+        outdoor_temperature = self.building.external_temperature[cte.HOUR][
+                              design_period_start_index:design_period_end_index]
+        results = HeatPumpBoilerTesHeating(hp=hp,
+                                           boiler=boiler,
+                                           tes=tes,
+                                           hourly_heating_demand_joules=heating_demand_joules,
+                                           heating_peak_load_watts=heating_peak_load_watts,
+                                           upper_limit_tes=upper_limit_tes_heating,
+                                           outdoor_temperature=outdoor_temperature,
+                                           dt=self.dt).simulation()
+        if min(results['TES Temperature']) < 35:
+          self.feasibility = False
+      elif cte.DOMESTIC_HOT_WATER in self.demand_types:
+        hp = self.energy_system.generation_systems[0]
+        hp.nominal_heat_output = self.individual['Generation Components'][0]['heating_capacity']
+        tes = self.energy_system.generation_systems[0].energy_storage_systems[0]
+        tes.volume = self.individual['Energy Storage Components'][0]['volume']
+        tes.height = self.building.average_storey_height - 1
+        tes.heating_coil_capacity = self.individual['Energy Storage Components'][0]['heating_coil_capacity'] \
+          if self.individual['Energy Storage Components'][0]['heating_coil_capacity'] is not None else None
+        dhw_demand_joules = self.design_period_energy_demands[cte.DOMESTIC_HOT_WATER]['demands']
+        upper_limit_tes = 65
+        design_period_start_index = self.design_period_energy_demands[cte.DOMESTIC_HOT_WATER]['start_index']
+        design_period_end_index = self.design_period_energy_demands[cte.DOMESTIC_HOT_WATER]['end_index']
+        outdoor_temperature = self.building.external_temperature[cte.HOUR][
+                              design_period_start_index:design_period_end_index]
+        results = DomesticHotWaterHeatPumpTes(hp=hp,
+                                              tes=tes,
+                                              hourly_dhw_demand_joules=dhw_demand_joules,
+                                              upper_limit_tes=upper_limit_tes,
+                                              outdoor_temperature=outdoor_temperature,
+                                              dt=self.dt).simulation()
+        if min(results['DHW TES Temperature']) < 55:
+          self.feasibility = False
+      if self.feasibility:
+        generation_system_types = [generation_system.system_type for generation_system in
+                                   self.energy_system.generation_systems]
+        for generation_component in self.individual['Generation Components']:
+          if generation_component['type'] in generation_system_types:
+            index = generation_system_types.index(generation_component['type'])
+            for demand_type in self.demand_types:
+              if demand_type in self.energy_system.generation_systems[index].energy_consumption:
+                generation_component['total_energy_consumption(kWh)'] = (sum(
+                  self.energy_system.generation_systems[index].energy_consumption[demand_type][cte.HOUR]) / 3.6e6)
+        for storage_component in self.individual['Energy Storage Components']:
+          if storage_component['type'] == f'{cte.THERMAL}_storage' and storage_component[
+            'heating_coil_capacity'] is not None:
+            for generation_system in self.energy_system.generation_systems:
+              if generation_system.energy_storage_systems is not None:
+                for storage_system in generation_system.energy_storage_systems:
+                  if storage_system.type_energy_stored == cte.THERMAL:
+                    for demand_type in self.demand_types:
+                      if demand_type in storage_system.heating_coil_energy_consumption:
+                        storage_component['total_energy_consumption(kWh)'] = (sum(
+                          storage_system.heating_coil_energy_consumption[demand_type][cte.HOUR]) / 3.6e6)
+
+  def life_cycle_cost_calculation(self, investment_cost, operation_cost_year_0, maintenance_cost_year_0,
+                                  life_cycle_duration=41):
+    """
+    Calculating the life cycle cost of the energy system. The original cost workflow in hub is not used to reduce
+    computation time.Here are the steps:
+        1- Costs catalog and different components are imported.
+        2- Capital costs (investment and reposition) are calculated and appended to a list to have the capital cash
+           flow.
+        3-
+    :param maintenance_cost_year_0:
+    :param operation_cost_year_0:
+    :param investment_cost:
+    :param life_cycle_duration:
+    :return:
+    """
+    capital_costs_cash_flow = [investment_cost]
+    operational_costs_cash_flow = [0]
+    maintenance_costs_cash_flow = [0]
+    end_of_life_costs = [0] * (life_cycle_duration + 1)
+    for i in range(1, life_cycle_duration + 1):
+      yearly_operational_cost = math.pow(1 + self.fuel_price_index, i) * operation_cost_year_0
+      yearly_maintenance_cost = math.pow(1 + self.consumer_price_index, i) * maintenance_cost_year_0
+      yearly_capital_cost = 0
+      for generation_system in self.individual['Generation Components']:
+        if (i % generation_system['lifetime']) == 0 and i != (life_cycle_duration - 1):
+          cost_increase = math.pow(1 + self.consumer_price_index, i)
+          heating_capacity = 0 if generation_system['heating_capacity'] is None else generation_system[
+            'heating_capacity']
+          cooling_capacity = 0 if generation_system['cooling_capacity'] is None else generation_system[
+            'cooling_capacity']
+          capacity = max(heating_capacity, cooling_capacity)
+          yearly_capital_cost += generation_system['unit_reposition_cost'] * capacity * cost_increase / 1000
+          yearly_capital_cost += -npf.pmt(self.interest_rate, self.credit_years,
+                                          investment_cost * self.percentage_credit)
+      for storage_system in self.individual['Energy Storage Components']:
+        if (i % storage_system['lifetime']) == 0 and i != (life_cycle_duration - 1):
+          cost_increase = math.pow(1 + self.consumer_price_index, i)
+          capacity = storage_system['volume'] if cte.THERMAL in storage_system['type'] else storage_system['capacity']
+          yearly_capital_cost += storage_system['unit_reposition_cost'] * capacity * cost_increase
+          yearly_capital_cost += -npf.pmt(self.interest_rate, self.credit_years,
+                                          investment_cost * self.percentage_credit)
+      capital_costs_cash_flow.append(yearly_capital_cost)
+      operational_costs_cash_flow.append(yearly_operational_cost)
+      maintenance_costs_cash_flow.append(yearly_maintenance_cost)
+    for year in range(1, life_cycle_duration + 1):
+      price_increase = math.pow(1 + self.consumer_price_index, year)
+      if year == life_cycle_duration:
+        end_of_life_costs[year] = (
+            self.building.thermal_zones_from_internal_zones[0].total_floor_area *
+            self.individual['End of Life Cost'] * price_increase
+        )
+
+    life_cycle_capital_cost = npf.npv(self.discount_rate, capital_costs_cash_flow)
+    life_cycle_operational_cost = npf.npv(self.discount_rate, operational_costs_cash_flow)
+    life_cycle_maintenance_cost = npf.npv(self.discount_rate, maintenance_costs_cash_flow)
+    life_cycle_end_of_life_cost = npf.npv(self.discount_rate, end_of_life_costs)
+    total_life_cycle_cost = life_cycle_capital_cost + life_cycle_operational_cost + life_cycle_maintenance_cost + life_cycle_end_of_life_cost
+    return total_life_cycle_cost
+
+  def costs_archetype(self):
+    costs_catalogue = CostsCatalogFactory('montreal_new').catalog.entries().archetypes
+    dictionary = Dictionaries().hub_function_to_montreal_custom_costs_function
+    costs_archetype = None
+    for archetype in costs_catalogue:
+      if dictionary[str(self.building.function)] == str(archetype.function):
+        costs_archetype = archetype
+    return costs_archetype
+
+  def unit_investment_cost(self, component_category, component_type):
+    """
+    Reading the investment and reposition costs of any component from costs catalogue
+    :param component_category: Due to the categorizations in the cost catalogue, we need this parameter to realize if
+    the component is a generation or storage component
+    :param component_type: Type of the component
+    :return:
+    """
+    investment_cost = 0
+    reposition_cost = 0
+    lifetime = 0
+    name = ''
+    capital_costs_chapter = self.costs.capital_cost.chapter('D_services')
+    if component_category == 'Generation':
+      generation_systems = self.energy_system.generation_systems
+      for generation_system in generation_systems:
+        if component_type == generation_system.system_type:
+          if generation_system.system_type == cte.HEAT_PUMP:
+            name += (
+                generation_system.source_medium.lower() + '_to_' + generation_system.supply_medium.lower() + '_' +
+                generation_system.system_type.lower().replace(' ', '_'))
+          elif generation_system.system_type == cte.BOILER:
+            if generation_system.fuel_type == cte.ELECTRICITY:
+              name += cte.ELECTRICAL.lower() + f'_{generation_system.system_type}'.lower()
+            else:
+              name += generation_system.fuel_type.lower() + f'_{generation_system.system_type}'.lower()
+          elif generation_system.system_type == cte.PHOTOVOLTAIC:
+            name += 'D2010_photovoltaic_system'
+          else:
+            if cte.HEATING or cte.DOMESTIC_HOT_WATER in self.demand_types:
+              name += 'template_heat'
+            else:
+              name += 'template_cooling'
+      for item in capital_costs_chapter.items:
+        if name in item.type:
+          investment_cost += float(capital_costs_chapter.item(item.type).initial_investment[0])
+          reposition_cost += float(capital_costs_chapter.item(item.type).reposition[0])
+          lifetime += float(capital_costs_chapter.item(item.type).lifetime)
+    elif component_category == 'Storage':
+      for generation_system in self.energy_system.generation_systems:
+        if generation_system.energy_storage_systems is not None:
+          for energy_storage_system in generation_system.energy_storage_systems:
+            if energy_storage_system.type_energy_stored == cte.THERMAL:
+              if energy_storage_system.heating_coil_capacity is not None:
+                investment_cost += float(capital_costs_chapter.item('D306010_storage_tank').initial_investment[0])
+                reposition_cost += float(capital_costs_chapter.item('D306010_storage_tank').reposition[0])
+                lifetime += float(capital_costs_chapter.item('D306010_storage_tank').lifetime)
+              else:
+                investment_cost += float(
+                  capital_costs_chapter.item('D306020_storage_tank_with_coil').initial_investment[0])
+                reposition_cost += float(capital_costs_chapter.item('D306020_storage_tank_with_coil').reposition[0])
+                lifetime += float(capital_costs_chapter.item('D306020_storage_tank_with_coil').lifetime)
+
+    return investment_cost, reposition_cost, lifetime
+
+  def unit_maintenance_cost(self, generation_system):
+    hvac_maintenance = self.costs.operational_cost.maintenance_hvac
+    pv_maintenance = self.costs.operational_cost.maintenance_pv
+    maintenance_cost = 0
+    component = None
+    if generation_system.system_type == cte.HEAT_PUMP and generation_system.source_medium == cte.AIR:
+      component = 'air_source_heat_pump'
+    elif generation_system.system_type == cte.HEAT_PUMP and generation_system.source_medium == cte.GROUND:
+      component = 'ground_source_heat_pump'
+    elif generation_system.system_type == cte.HEAT_PUMP and generation_system.source_medium == cte.WATER:
+      component = 'water_source_heat_pump'
+    elif generation_system.system_type == cte.BOILER and generation_system.fuel_type == cte.GAS:
+      component = 'gas_boiler'
+    elif generation_system.system_type == cte.BOILER and generation_system.fuel_type == cte.ELECTRICITY:
+      component = 'electric_boiler'
+    elif generation_system.system_type == cte.PHOTOVOLTAIC:
+      maintenance_cost += pv_maintenance
+    else:
+      if cte.HEATING or cte.DOMESTIC_HOT_WATER in self.demand_types:
+        component = 'general_heating_equipment'
+      else:
+        component = 'general_cooling_equipment'
+    for item in hvac_maintenance:
+      if item.type == component:
+        maintenance_cost += item.maintenance[0]
+    return maintenance_cost
+
+  def fuel_cost_per_kwh(self, fuel_type, fuel_cost_tariff_type):
+    fuel_cost = 0
+    catalog_fuels = self.costs.operational_cost.fuels
+    for fuel in catalog_fuels:
+      if fuel_type == fuel.type and fuel_cost_tariff_type == fuel.variable.rate_type:
+        if fuel.type == cte.ELECTRICITY and fuel_cost_tariff_type == 'fixed':
+          fuel_cost = fuel.variable.values[0]
+        elif fuel.type == cte.ELECTRICITY and fuel_cost_tariff_type == 'variable':
+          fuel_cost = fuel.variable.values[0]
+        else:
+          if fuel.type == cte.BIOMASS:
+            conversion_factor = 1
+          else:
+            conversion_factor = fuel.density[0]
+          fuel_cost = fuel.variable.values[0] / (conversion_factor * fuel.lower_heating_value[0] * 0.277)
+    return fuel_cost

energy_system_modelling_package/energy_system_modelling_factories/system_sizing_methods/genetic_algorithm/multi_objective_genetic_algorithm.py

@@ -0,0 +1,418 @@
+import copy
+import math
+import random
+import hub.helpers.constants as cte
+from energy_system_modelling_package.energy_system_modelling_factories.system_sizing_methods.genetic_algorithm.individual import \
+  Individual
+import matplotlib.pyplot as plt
+
+
+class MultiObjectiveGeneticAlgorithm:
+  def __init__(self, population_size=100, generations=50, crossover_rate=0.8, mutation_rate=0.1,
+               optimization_scenario=None, output_path=None):
+    self.population_size = population_size
+    self.population = []
+    self.generations = generations
+    self.crossover_rate = crossover_rate
+    self.mutation_rate = mutation_rate
+    self.optimization_scenario = optimization_scenario
+    self.list_of_solutions = []
+    self.best_solution = None
+    self.best_solution_generation = None
+    self.output_path = output_path
+
+# Initialize Population
+  def initialize_population(self, building, energy_system):
+    design_period_energy_demands = self.design_period_identification(building)
+    attempts = 0
+    max_attempts = self.population_size * 5
+    while len(self.population) < self.population_size and attempts < max_attempts:
+      individual = Individual(building=building,
+                              energy_system=energy_system,
+                              design_period_energy_demands=design_period_energy_demands,
+                              optimization_scenario=self.optimization_scenario)
+      individual.initialization()
+      attempts += 1
+      if self.initial_population_feasibility_check(individual, energy_system.demand_types,
+                                                   design_period_energy_demands):
+        self.population.append(individual)
+    if len(self.population) < self.population_size:
+      raise RuntimeError(f"Could not generate a feasible population of size {self.population_size}. "
+                         f"Only {len(self.population)} feasible individuals were generated.")
+
+  @staticmethod
+  def initial_population_feasibility_check(individual, demand_types, design_period_demands):
+    total_heating_capacity = sum(
+      component['heating_capacity'] for component in individual.individual['Generation Components']
+      if component['heating_capacity'] is not None
+    )
+    total_cooling_capacity = sum(
+      component['cooling_capacity'] for component in individual.individual['Generation Components']
+      if component['cooling_capacity'] is not None
+    )
+    max_heating_demand = max(design_period_demands[cte.HEATING]['demands']) / cte.WATTS_HOUR_TO_JULES
+    max_cooling_demand = max(design_period_demands[cte.COOLING]['demands']) / cte.WATTS_HOUR_TO_JULES
+    max_dhw_demand = max(design_period_demands[cte.DOMESTIC_HOT_WATER]['demands']) / cte.WATTS_HOUR_TO_JULES
+    if cte.HEATING in demand_types and total_heating_capacity < 0.5 * max_heating_demand:
+      return False
+    if cte.DOMESTIC_HOT_WATER in demand_types and total_heating_capacity < 0.5 * max_dhw_demand:
+      return False
+    if cte.COOLING in demand_types and total_cooling_capacity < 0.5 * max_cooling_demand:
+      return False
+    total_volume = sum(
+      component['volume'] for component in individual.individual['Energy Storage Components']
+      if component['volume'] is not None
+    )
+    max_storage_volume = individual.available_space * 0.1
+    if total_volume > max_storage_volume:
+      return False
+    return True
+
+  def nsga2_selection(self, population, fronts, crowding_distances):
+    new_population = []
+    i = 0
+    while len(new_population) + len(fronts[i]) <= self.population_size:
+      for index in fronts[i]:
+        # Skip individuals with infinite fitness values to avoid propagation
+        if not math.isinf(self.population[index].individual['fitness_score'][0]) and \
+            not math.isinf(self.population[index].individual['fitness_score'][1]):
+          new_population.append(population[index])
+      i += 1
+      if i >= len(fronts):
+        break
+    if len(new_population) < self.population_size and i < len(fronts):
+      sorted_front = sorted(fronts[i], key=lambda x: crowding_distances[i][x], reverse=True)
+      for index in sorted_front:
+        if len(new_population) < self.population_size:
+          if not math.isinf(self.population[index].individual['fitness_score'][0]) and \
+              not math.isinf(self.population[index].individual['fitness_score'][1]):
+            new_population.append(population[index])
+        else:
+          break
+    return new_population
+
+  def fast_non_dominated_sort(self):
+    population_size = self.population_size
+    s = [[] for _ in range(population_size)]
+    front = [[]]
+    n = [0] * population_size
+    rank = [0] * population_size
+    for p in range(population_size):
+      s[p] = []
+      n[p] = 0
+      for q in range(population_size):
+        if self.dominates(self.population[p], self.population[q]):
+          s[p].append(q)
+        elif self.dominates(self.population[q], self.population[p]):
+          n[p] += 1
+      if n[p] == 0:
+        rank[p] = 0
+        front[0].append(p)
+    i = 0
+    while front[i]:
+      next_front = set()
+      for p in front[i]:
+        for q in s[p]:
+          n[q] -= 1
+          if n[q] == 0:
+            rank[q] = i + 1
+            next_front.add(q)
+      i += 1
+      front.append(list(next_front))
+    del front[-1]
+    return front
+
+  @staticmethod
+  def dominates(individual1, individual2):
+    lcc1, lce1 = individual1.individual['fitness_score']
+    lcc2, lce2 = individual2.individual['fitness_score']
+    return (lcc1 <= lcc2 and lce1 <= lce2) and (lcc1 < lcc2 or lce1 < lce2)
+
+  def calculate_crowding_distance(self, front):
+    crowding_distance = [0] * len(self.population)
+    for objective in ['lcc', 'total_energy_consumption']:
+      sorted_front = sorted(front, key=lambda x: self.population[x].individual[objective])
+      # Set distances to finite large numbers rather than `inf`
+      crowding_distance[sorted_front[0]] = float(1e9)
+      crowding_distance[sorted_front[-1]] = float(1e9)
+      objective_min = self.population[sorted_front[0]].individual[objective]
+      objective_max = self.population[sorted_front[-1]].individual[objective]
+      if objective_max != objective_min:
+        for i in range(1, len(sorted_front) - 1):
+          crowding_distance[sorted_front[i]] += (
+              (self.population[sorted_front[i + 1]].individual[objective] -
+               self.population[sorted_front[i - 1]].individual[objective]) /
+              (objective_max - objective_min))
+    return crowding_distance
+
+  def crossover(self, parent1, parent2):
+    """
+    Crossover between two parents to produce two children.
+    swaps generation components and storage components between the two parents with a 50% chance.
+    :param parent1: First parent individual.
+    :param parent2: second parent individual.
+    :return: Two child individuals (child1 and child2).
+    """
+    if random.random() < self.crossover_rate:
+      # Deep copy of the parents to create children
+      child1, child2 = copy.deepcopy(parent1), copy.deepcopy(parent2)
+      # Crossover for Generation Components
+      for i in range(len(parent1.individual['Generation Components'])):
+        if random.random() < 0.5:
+          # swap the entire generation component
+          child1.individual['Generation Components'][i], child2.individual['Generation Components'][i] = (
+            child2.individual['Generation Components'][i],
+            child1.individual['Generation Components'][i]
+          )
+
+      # Crossover for Energy storage Components
+      for i in range(len(parent1.individual['Energy Storage Components'])):
+        if random.random() < 0.5:
+          # swap the entire storage component
+          child1.individual['Energy Storage Components'][i], child2.individual['Energy Storage Components'][i] = (
+            child2.individual['Energy Storage Components'][i],
+            child1.individual['Energy Storage Components'][i]
+          )
+
+      return child1, child2
+    else:
+      # If crossover does not happen, return copies of the original parents
+      return copy.deepcopy(parent1), copy.deepcopy(parent2)
+
+  def mutate(self, individual, building, energy_system):
+    """
+    Mutates the individual's generation and storage components.
+
+    - `individual`: The individual to mutate (contains generation and storage components).
+    - `building`: Building data that contains constraints such as peak heating load and available space.
+
+    Returns the mutated individual.
+    """
+    design_period_energy_demands = self.design_period_identification(building)
+    # Mutate Generation Components
+    for generation_component in individual['Generation Components']:
+      if random.random() < self.mutation_rate:
+        if (generation_component['nominal_heating_efficiency'] is not None and cte.HEATING or cte.DOMESTIC_HOT_WATER in
+            energy_system.demand_types):
+          # Mutate heating capacity
+          if cte.HEATING in energy_system.demand_types:
+            generation_component['heating_capacity'] = random.uniform(
+              0, max(design_period_energy_demands[cte.HEATING]['demands']) / cte.WATTS_HOUR_TO_JULES)
+          else:
+            generation_component['heating_capacity'] = random.uniform(
+              0, max(design_period_energy_demands[cte.DOMESTIC_HOT_WATER]['demands']) / cte.WATTS_HOUR_TO_JULES)
+        if generation_component['nominal_cooling_efficiency'] is not None and cte.COOLING in energy_system.demand_types:
+          # Mutate cooling capacity
+          generation_component['cooling_capacity'] = random.uniform(
+            0, max(design_period_energy_demands[cte.COOLING]['demands']) / cte.WATTS_HOUR_TO_JULES)
+    # Mutate storage Components
+    for storage_component in individual['Energy Storage Components']:
+      if random.random() < self.mutation_rate:
+        if storage_component['type'] == f'{cte.THERMAL}_storage':
+          # Mutate the volume of thermal storage
+          max_available_space = 0.01 * building.volume / building.storeys_above_ground
+          storage_component['volume'] = random.uniform(0, max_available_space)
+          if storage_component['heating_coil_capacity'] is not None:
+            if cte.HEATING in energy_system.demand_types:
+              storage_component['heating_coil_capacity'] = random.uniform(0, max(
+                design_period_energy_demands[cte.HEATING]['demands']) / cte.WATTS_HOUR_TO_JULES)
+            else:
+              storage_component['heating_coil_capacity'] = random.uniform(0, max(
+                design_period_energy_demands[cte.DOMESTIC_HOT_WATER]['demands']) / cte.WATTS_HOUR_TO_JULES)
+    return individual
+
+  def solve_ga(self, building, energy_system):
+    self.initialize_population(building, energy_system)
+    for individual in self.population:
+      individual.score_evaluation()
+    pareto_population = []
+
+    for generation in range(1, self.generations + 1):
+      print(f"Generation {generation}")
+      fronts = self.fast_non_dominated_sort()
+
+      # Ensure the front calculation is valid
+      if not fronts or not fronts[0]:
+        print("Warning: No valid non-dominated front found.")
+        continue
+
+      # Calculate crowding distances and select individuals
+      crowding_distances = [self.calculate_crowding_distance(front) for front in fronts]
+      self.population = self.nsga2_selection(self.population, fronts, crowding_distances)
+
+      # Add only valid indices to pareto_population
+      pareto_population.extend([self.population[i] for i in fronts[0] if i < len(self.population)])
+
+      # Check population bounds
+      if len(pareto_population) > self.population_size:
+        pareto_population = pareto_population[:self.population_size]
+
+      # Generate the next population through crossover and mutation
+      next_population = []
+      while len(next_population) < self.population_size:
+        parent1, parent2 = random.choice(self.population), random.choice(self.population)
+        child1, child2 = self.crossover(parent1, parent2)
+        self.mutate(child1.individual, building, energy_system)
+        self.mutate(child2.individual, building, energy_system)
+        child1.score_evaluation()
+        child2.score_evaluation()
+        next_population.extend([child1, child2])
+
+      self.population = next_population[:self.population_size]
+
+    # Ensure pareto_population contains the correct non-dominated individuals before TOPSIS
+    if not pareto_population:
+      print("No Pareto solutions found during optimization.")
+      return None
+
+    # Recalculate pareto front indices based on updated pareto_population
+    fronts = self.fast_non_dominated_sort()
+    pareto_front_indices = fronts[0] if fronts else []
+    pareto_front_indices = [i for i in pareto_front_indices if i < len(pareto_population)]
+    print(f"Final pareto_front_indices: {pareto_front_indices}, pareto_population size: {len(pareto_population)}")
+
+    if not pareto_front_indices:
+      print("No valid solution found after TOPSIS due to empty pareto front indices.")
+      return None
+
+    global_pareto_front = [pareto_population[i] for i in pareto_front_indices]
+
+    # Get the top N solutions with TOPSIS
+    top_n = 3  # Adjust this value based on how many top solutions you want to explore
+    self.best_solution = self.topsis_decision_making(global_pareto_front, top_n=top_n)
+
+    # Print the top N solutions
+    if self.best_solution:
+      print("Top solutions after TOPSIS:")
+      for idx, solution in enumerate(self.best_solution, 1):
+        print(f"Solution {idx}: LCC = {solution.individual['lcc']}, "
+              f"LCE = {solution.individual['total_energy_consumption']}")
+    else:
+      print("No valid solutions found after TOPSIS.")
+
+    if pareto_population:
+      self.plot_pareto_front(pareto_population)
+
+    return self.best_solution
+
+  @staticmethod
+  def design_period_identification(building):
+    def get_start_end_indices(max_day_index, total_days):
+      if max_day_index > 0 and max_day_index < total_days - 1:
+        start_index = (max_day_index - 1) * 24
+        end_index = (max_day_index + 2) * 24
+      elif max_day_index == 0:
+        start_index = 0
+        end_index = (max_day_index + 2) * 24
+      else:
+        start_index = (max_day_index - 1) * 24
+        end_index = total_days * 24
+      return start_index, end_index
+
+    # Calculate daily demands
+    heating_daily_demands = [sum(building.heating_demand[cte.HOUR][i:i + 24]) for i in
+                             range(0, len(building.heating_demand[cte.HOUR]), 24)]
+    cooling_daily_demands = [sum(building.cooling_demand[cte.HOUR][i:i + 24]) for i in
+                             range(0, len(building.cooling_demand[cte.HOUR]), 24)]
+    dhw_daily_demands = [sum(building.domestic_hot_water_heat_demand[cte.HOUR][i:i + 24]) for i in
+                         range(0, len(building.domestic_hot_water_heat_demand[cte.HOUR]), 24)]
+    # Get the day with maximum demand for each type
+    heating_max_day = heating_daily_demands.index(max(heating_daily_demands))
+    cooling_max_day = cooling_daily_demands.index(max(cooling_daily_demands))
+    dhw_max_day = dhw_daily_demands.index(max(dhw_daily_demands))
+    # Get the start and end indices for each demand type
+    heating_start, heating_end = get_start_end_indices(heating_max_day, len(heating_daily_demands))
+    cooling_start, cooling_end = get_start_end_indices(cooling_max_day, len(cooling_daily_demands))
+    dhw_start, dhw_end = get_start_end_indices(dhw_max_day, len(dhw_daily_demands))
+    # Return the design period energy demands
+    return {
+      f'{cte.HEATING}': {'demands': building.heating_demand[cte.HOUR][heating_start:heating_end],
+                         'start_index': heating_start, 'end_index': heating_end},
+      f'{cte.COOLING}': {'demands': building.cooling_demand[cte.HOUR][cooling_start:cooling_end],
+                         'start_index': cooling_start, 'end_index': cooling_end},
+      f'{cte.DOMESTIC_HOT_WATER}': {'demands': building.domestic_hot_water_heat_demand[cte.HOUR][dhw_start:dhw_end],
+                                    'start_index': dhw_start, 'end_index': dhw_end}
+    }
+
+  def topsis_decision_making(self, pareto_front, top_n=5):
+    """
+    Perform TOPSIS decision-making to choose the best solutions from the Pareto front.
+
+    :param pareto_front: List of individuals in the Pareto front
+    :param top_n: Number of top solutions to select based on TOPSIS ranking
+    :return: List of top N individuals based on TOPSIS ranking
+    """
+    # Filter out infinite values from the pareto front before processing
+    pareto_front = [ind for ind in pareto_front if all(math.isfinite(v) for v in ind.individual['fitness_score'])]
+    if not pareto_front:
+      return None
+
+    # Step 1: Normalize the objective functions (cost and energy consumption)
+    min_lcc = min([ind.individual['lcc'] for ind in pareto_front])
+    max_lcc = max([ind.individual['lcc'] for ind in pareto_front])
+    min_lce = min([ind.individual['total_energy_consumption'] for ind in pareto_front])
+    max_lce = max([ind.individual['total_energy_consumption'] for ind in pareto_front])
+
+    normalized_pareto_front = []
+    for ind in pareto_front:
+      normalized_lcc = (ind.individual['lcc'] - min_lcc) / (max_lcc - min_lcc) if max_lcc > min_lcc else 0
+      normalized_lce = (ind.individual['total_energy_consumption'] - min_lce) / (
+          max_lce - min_lce) if max_lce > min_lce else 0
+      normalized_pareto_front.append((ind, normalized_lcc, normalized_lce))
+
+    # Step 2: Calculate the ideal and worst solutions
+    ideal_solution = (0, 0)  # Ideal is minimum LCC and minimum LCE (0, 0 after normalization)
+    worst_solution = (1, 1)  # Worst is maximum LCC and maximum LCE (1, 1 after normalization)
+
+    # Step 3: Calculate the distance to the ideal and worst solutions
+    best_distances = []
+    worst_distances = []
+    for ind, normalized_lcc, normalized_lce in normalized_pareto_front:
+      distance_to_ideal = math.sqrt(
+        (normalized_lcc - ideal_solution[0]) ** 2 + (normalized_lce - ideal_solution[1]) ** 2)
+      distance_to_worst = math.sqrt(
+        (normalized_lcc - worst_solution[0]) ** 2 + (normalized_lce - worst_solution[1]) ** 2)
+      best_distances.append(distance_to_ideal)
+      worst_distances.append(distance_to_worst)
+
+    # Step 4: Calculate relative closeness to the ideal solution
+    similarity = [worst / (best + worst) for best, worst in zip(best_distances, worst_distances)]
+
+    # Step 5: Select the top N individuals with the highest similarity scores
+    top_indices = sorted(range(len(similarity)), key=lambda i: similarity[i], reverse=True)[:top_n]
+    top_solutions = [pareto_front[i] for i in top_indices]
+
+    # Plot the similarity scores
+    self.plot_topsis_similarity(similarity)
+
+    return top_solutions
+
+  @staticmethod
+  def plot_topsis_similarity(similarity):
+    """
+    Plot the TOPSIS similarity scores for visualization.
+
+    :param similarity: List of similarity scores for each individual in the Pareto front
+    """
+    plt.figure(figsize=(10, 6))
+    plt.plot(range(len(similarity)), similarity, 'bo-', label='TOPSIS Similarity Scores')
+    plt.xlabel('Pareto Front Solution Index')
+    plt.ylabel('Similarity Score')
+    plt.title('TOPSIS Similarity Scores for Pareto Front Solutions')
+    plt.legend()
+    plt.grid(True)
+    plt.show()
+
+  @staticmethod
+  def plot_pareto_front(pareto_population):
+    # Extract LCC and LCE for plotting
+    lcc_values = [individual.individual['lcc'] for individual in pareto_population]
+    lce_values = [individual.individual['total_energy_consumption'] for individual in pareto_population]
+    plt.figure(figsize=(10, 6))
+    plt.scatter(lcc_values, lce_values, color='blue', label='Pareto Front', alpha=0.6, edgecolors='w', s=80)
+    plt.title('Pareto Front for Life Cycle Cost vs Life Cycle Energy')
+    plt.xlabel('Life Cycle Cost (LCC)')
+    plt.ylabel('Life Cycle Energy (LCE)')
+    plt.grid(True)
+    plt.legend()
+    plt.show()

energy_system_modelling_package/energy_system_modelling_factories/system_sizing_methods/genetic_algorithm/multi_objective_genetic_algorithm_rethinking.py

@@ -0,0 +1,94 @@
+import copy
+import math
+import random
+import hub.helpers.constants as cte
+from energy_system_modelling_package.energy_system_modelling_factories.system_sizing_methods.genetic_algorithm.individual import \
+  Individual
+import matplotlib.pyplot as plt
+import time
+
+
+class MultiObjectiveGeneticAlgorithm:
+  def __init__(self, population_size=100, generations=50, crossover_rate=0.8, mutation_rate=0.1,
+               optimization_scenario=None, output_path=None):
+    self.population_size = population_size
+    self.population = []
+    self.generations = generations
+    self.crossover_rate = crossover_rate
+    self.mutation_rate = mutation_rate
+    self.optimization_scenario = optimization_scenario
+    self.list_of_solutions = []
+    self.best_solution = None
+    self.best_solution_generation = None
+    self.output_path = output_path
+
+# Initialize Population
+  def initialize_population(self, building, energy_system):
+    design_period_start_time = time.time()
+    design_period_energy_demands = self.design_period_identification(building)
+    design_period_time = time.time() - design_period_start_time
+    print(f"design period identification took {design_period_time:.2f} seconds")
+    initializing_time_start = time.time()
+    attempts = 0
+    max_attempts = self.population_size * 20
+    while len(self.population) < self.population_size and attempts < max_attempts:
+      individual = Individual(building=building,
+                              energy_system=energy_system,
+                              design_period_energy_demands=design_period_energy_demands,
+                              optimization_scenario=self.optimization_scenario)
+      individual.initialization()
+      individual.score_evaluation()
+      attempts += 1
+      if individual.feasibility:
+        self.population.append(individual)
+    if len(self.population) < self.population_size:
+      raise RuntimeError(f"Could not generate a feasible population of size {self.population_size}. "
+                         f"Only {len(self.population)} feasible individuals were generated.")
+    initializing_time = time.time() - initializing_time_start
+    print(f"initializing took {initializing_time:.2f} seconds")
+
+
+  @staticmethod
+  def design_period_identification(building):
+    def get_start_end_indices(max_day_index, total_days):
+      if max_day_index > 0 and max_day_index < total_days - 1:
+        start_index = (max_day_index - 1) * 24
+        end_index = (max_day_index + 2) * 24
+      elif max_day_index == 0:
+        start_index = 0
+        end_index = (max_day_index + 2) * 24
+      else:
+        start_index = (max_day_index - 1) * 24
+        end_index = total_days * 24
+      return start_index, end_index
+
+    # Calculate daily demands
+    heating_daily_demands = [sum(building.heating_demand[cte.HOUR][i:i + 24]) for i in
+                             range(0, len(building.heating_demand[cte.HOUR]), 24)]
+    cooling_daily_demands = [sum(building.cooling_demand[cte.HOUR][i:i + 24]) for i in
+                             range(0, len(building.cooling_demand[cte.HOUR]), 24)]
+    dhw_daily_demands = [sum(building.domestic_hot_water_heat_demand[cte.HOUR][i:i + 24]) for i in
+                         range(0, len(building.domestic_hot_water_heat_demand[cte.HOUR]), 24)]
+    # Get the day with maximum demand for each type
+    heating_max_day = heating_daily_demands.index(max(heating_daily_demands))
+    cooling_max_day = cooling_daily_demands.index(max(cooling_daily_demands))
+    dhw_max_day = dhw_daily_demands.index(max(dhw_daily_demands))
+    # Get the start and end indices for each demand type
+    heating_start, heating_end = get_start_end_indices(heating_max_day, len(heating_daily_demands))
+    cooling_start, cooling_end = get_start_end_indices(cooling_max_day, len(cooling_daily_demands))
+    dhw_start, dhw_end = get_start_end_indices(dhw_max_day, len(dhw_daily_demands))
+    # Return the design period energy demands
+    return {
+      f'{cte.HEATING}': {'demands': building.heating_demand[cte.HOUR][heating_start:heating_end],
+                         'start_index': heating_start, 'end_index': heating_end},
+      f'{cte.COOLING}': {'demands': building.cooling_demand[cte.HOUR][cooling_start:cooling_end],
+                         'start_index': cooling_start, 'end_index': cooling_end},
+      f'{cte.DOMESTIC_HOT_WATER}': {'demands': building.domestic_hot_water_heat_demand[cte.HOUR][dhw_start:dhw_end],
+                                    'start_index': dhw_start, 'end_index': dhw_end}
+    }
+
+  def solve_ga(self, building, energy_system):
+    self.initialize_population(building, energy_system)
+    for individual in self.population:
+      print(individual.individual)
+    print([ind.fitness_score for ind in self.population])

energy_system_modelling_package/energy_system_modelling_factories/system_sizing_methods/genetic_algorithm/single_objective_genetic_algorithm.py

@@ -0,0 +1,342 @@
+import copy
+import math
+import random
+import hub.helpers.constants as cte
+from energy_system_modelling_package.energy_system_modelling_factories.system_sizing_methods.genetic_algorithm.individual import \
+  Individual
+
+
+class SingleObjectiveGeneticAlgorithm:
+  def __init__(self, population_size=100, generations=20, crossover_rate=0.8, mutation_rate=0.1,
+               optimization_scenario=None, output_path=None):
+    self.population_size = population_size
+    self.population = []
+    self.generations = generations
+    self.crossover_rate = crossover_rate
+    self.mutation_rate = mutation_rate
+    self.optimization_scenario = optimization_scenario
+    self.list_of_solutions = []
+    self.best_solution = None
+    self.best_solution_generation = None
+    self.output_path = output_path
+
+  # Initialize Population
+  def initialize_population(self, building, energy_system):
+    """
+    Initialize a population of individuals with feasible configurations for optimizing the sizes of
+    generation and storage components of an energy system.
+
+    :param building: Building object with associated data
+    :param energy_system: Energy system to optimize
+    """
+    design_period_energy_demands = self.design_period_identification(building)
+    attempts = 0  # Track attempts to avoid an infinite loop in rare cases
+    max_attempts = self.population_size * 5
+
+    while len(self.population) < self.population_size and attempts < max_attempts:
+      individual = Individual(building=building,
+                              energy_system=energy_system,
+                              design_period_energy_demands=design_period_energy_demands,
+                              optimization_scenario=self.optimization_scenario)
+
+      individual.initialization()
+      attempts += 1
+
+    # Enhanced feasibility check
+      if self.initial_population_feasibility_check(individual, energy_system.demand_types, design_period_energy_demands):
+        self.population.append(individual)
+
+    # Raise an error or print a warning if the population size goal is not met after max_attempts
+    if len(self.population) < self.population_size:
+      raise RuntimeError(f"Could not generate a feasible population of size {self.population_size}. "
+                         f"Only {len(self.population)} feasible individuals were generated.")
+
+  @staticmethod
+  def initial_population_feasibility_check(individual, demand_types, design_period_demands):
+    """
+    Check if the individual meets basic feasibility requirements for heating, cooling, and DHW capacities
+    and storage volume.
+
+    :param individual: Individual to check
+    :param demand_types: List of demand types (e.g., heating, cooling, DHW)
+    :param design_period_demands: Design period demand values for heating, cooling, and DHW
+    :return: True if feasible, False otherwise
+    """
+    # Calculate total heating and cooling capacities
+    total_heating_capacity = sum(
+      component['heating_capacity'] for component in individual.individual['Generation Components']
+      if component['heating_capacity'] is not None
+    )
+    total_cooling_capacity = sum(
+      component['cooling_capacity'] for component in individual.individual['Generation Components']
+      if component['cooling_capacity'] is not None
+    )
+
+    # Maximum demands for each demand type (converted to kW)
+    max_heating_demand = max(design_period_demands[cte.HEATING]['demands']) / cte.WATTS_HOUR_TO_JULES
+    max_cooling_demand = max(design_period_demands[cte.COOLING]['demands']) / cte.WATTS_HOUR_TO_JULES
+    max_dhw_demand = max(design_period_demands[cte.DOMESTIC_HOT_WATER]['demands']) / cte.WATTS_HOUR_TO_JULES
+
+    # Check heating capacity feasibility
+    if cte.HEATING in demand_types and total_heating_capacity < 0.5 * max_heating_demand:
+      return False
+
+    # Check DHW capacity feasibility
+    if cte.DOMESTIC_HOT_WATER in demand_types and total_heating_capacity < 0.5 * max_dhw_demand:
+      return False
+
+    # Check cooling capacity feasibility
+    if cte.COOLING in demand_types and total_cooling_capacity < 0.5 * max_cooling_demand:
+      return False
+
+    # Check storage volume feasibility
+    total_volume = sum(
+      component['volume'] for component in individual.individual['Energy Storage Components']
+      if component['volume'] is not None
+    )
+    # Limit storage to 10% of building's available space
+    max_storage_volume = individual.available_space * 0.1
+    if total_volume > max_storage_volume:
+      return False
+
+    return True  # Feasible if all checks are passed
+
+  def order_population(self):
+    """
+    ordering the population based on the fitness score in ascending order
+    :return:
+    """
+    self.population = sorted(self.population, key=lambda x: x.fitness_score)
+
+  def tournament_selection(self):
+    selected = []
+    for _ in range(len(self.population)):
+      i, j = random.sample(range(self.population_size), 2)
+      if self.population[i].individual['fitness_score'] < self.population[j].individual['fitness_score']:
+        selected.append(copy.deepcopy(self.population[i]))
+      else:
+        selected.append(copy.deepcopy(self.population[j]))
+    return selected
+
+  def crossover(self, parent1, parent2):
+    """
+    Crossover between two parents to produce two children.
+
+    swaps generation components and storage components between the two parents with a 50% chance.
+
+    :param parent1: First parent individual.
+    :param parent2: second parent individual.
+    :return: Two child individuals (child1 and child2).
+    """
+    if random.random() < self.crossover_rate:
+        # Deep copy of the parents to create children
+      child1, child2 = copy.deepcopy(parent1), copy.deepcopy(parent2)
+    # Crossover for Generation Components
+      for i in range(len(parent1.individual['Generation Components'])):
+        if random.random() < 0.5:
+            # swap the entire generation component
+          child1.individual['Generation Components'][i], child2.individual['Generation Components'][i] = (
+            child2.individual['Generation Components'][i],
+            child1.individual['Generation Components'][i]
+          )
+
+    # Crossover for Energy storage Components
+      for i in range(len(parent1.individual['Energy Storage Components'])):
+        if random.random() < 0.5:
+            # swap the entire storage component
+          child1.individual['Energy Storage Components'][i], child2.individual['Energy Storage Components'][i] = (
+            child2.individual['Energy Storage Components'][i],
+            child1.individual['Energy Storage Components'][i]
+          )
+
+      return child1, child2
+    else:
+        # If crossover does not happen, return copies of the original parents
+      return copy.deepcopy(parent1), copy.deepcopy(parent2)
+
+  def mutate(self, individual, building, energy_system):
+    """
+    Mutates the individual's generation and storage components.
+
+    - `individual`: The individual to mutate (contains generation and storage components).
+    - `building`: Building data that contains constraints such as peak heating load and available space.
+
+    Returns the mutated individual.
+    """
+    design_period_energy_demands = self.design_period_identification(building)
+    # Mutate Generation Components
+    for generation_component in individual['Generation Components']:
+      if random.random() < self.mutation_rate:
+        if (generation_component['nominal_heating_efficiency'] is not None and cte.HEATING or cte.DOMESTIC_HOT_WATER in
+                energy_system.demand_types):
+            # Mutate heating capacity
+          if cte.HEATING in energy_system.demand_types:
+            generation_component['heating_capacity'] = random.uniform(
+              0, max(design_period_energy_demands[cte.HEATING]['demands']) / cte.WATTS_HOUR_TO_JULES)
+          else:
+            generation_component['heating_capacity'] = random.uniform(
+              0, max(design_period_energy_demands[cte.DOMESTIC_HOT_WATER]['demands']) / cte.WATTS_HOUR_TO_JULES)
+        if generation_component['nominal_cooling_efficiency'] is not None and cte.COOLING in energy_system.demand_types:
+            # Mutate cooling capacity
+          generation_component['cooling_capacity'] = random.uniform(
+            0, max(design_period_energy_demands[cte.COOLING]['demands']) / cte.WATTS_HOUR_TO_JULES)
+    # Mutate storage Components
+    for storage_component in individual['Energy Storage Components']:
+      if random.random() < self.mutation_rate:
+        if storage_component['type'] == f'{cte.THERMAL}_storage':
+            # Mutate the volume of thermal storage
+          max_available_space = 0.01 * building.volume / building.storeys_above_ground
+          storage_component['volume'] = random.uniform(0, max_available_space)
+          if storage_component['heating_coil_capacity'] is not None:
+            if cte.HEATING in energy_system.demand_types:
+              storage_component['heating_coil_capacity'] = random.uniform(0, max(
+                design_period_energy_demands[cte.HEATING]['demands']) / cte.WATTS_HOUR_TO_JULES)
+            else:
+              storage_component['heating_coil_capacity'] = random.uniform(0, max(
+                design_period_energy_demands[cte.DOMESTIC_HOT_WATER]['demands']) / cte.WATTS_HOUR_TO_JULES)
+    return individual
+
+  def solve_ga(self, building, energy_system):
+    """
+    solving GA for a single energy system. Here are the steps:
+    1. Initialize population using the "initialize_population" method in this class.
+    2. Evaluate the initial population using the "score_evaluation" method in the Individual class.
+    3. sort population based on fitness score.
+    4. Repeat selection, crossover, and mutation for a fixed number of generations.
+    5. Track the best solution found during the optimization process.
+
+    :param building: Building object for the energy system.
+    :param energy_system: Energy system to optimize.
+    :return: Best solution after running the GA.
+    """
+    # step 1: Initialize the population
+    self.initialize_population(building, energy_system)
+    # step 2: Evaluate the initial population
+    for individual in self.population:
+      individual.score_evaluation()
+    # step 3: Order population based on fitness scores
+    self.order_population()
+    print([individual.fitness_score for individual in self.population])
+    # Track the best solution
+    self.best_solution = self.population[0]
+    self.best_solution_generation = 0
+    self.list_of_solutions.append(copy.deepcopy(self.best_solution.individual))
+    # step 4: Run GA for a fixed number of generations
+    for generation in range(1, self.generations):
+      print(f"Generation {generation}")
+    # selection (using tournament selection)
+      selected_population = self.tournament_selection()
+    # Create the next generation through crossover and mutation
+      next_population = []
+      for i in range(0, self.population_size, 2):
+        parent1 = selected_population[i]
+        parent2 = selected_population[i + 1] if (i + 1) < len(selected_population) else selected_population[0]
+    # step 5: Apply crossover
+        child1, child2 = self.crossover(parent1, parent2)
+    # step 6: Apply mutation
+        self.mutate(child1.individual, building, energy_system)
+        self.mutate(child2.individual, building, energy_system)
+    # step 7: Evaluate the children
+        child1.score_evaluation()
+        child2.score_evaluation()
+        next_population.extend([child1, child2])
+    # Replace old population with the new one
+      self.population = next_population
+    # step 8: sort the new population based on fitness
+      self.order_population()
+      print([individual.fitness_score for individual in self.population])
+    # Track the best solution found in this generation
+      if self.population[0].individual['fitness_score'] < self.best_solution.individual['fitness_score']:
+        self.best_solution = self.population[0]
+        self.best_solution_generation = generation
+    # store the best solution in the list of solutions
+      self.list_of_solutions.append(copy.deepcopy(self.population[0].individual))
+    print(f"Best solution found in generation {self.best_solution_generation}")
+    print(f"Best solution: {self.best_solution.individual}")
+    return self.best_solution
+
+  @staticmethod
+  def topsis_decision_making(pareto_front):
+    """
+    Perform TOPSIS decision-making to choose the best solution from the Pareto front.
+
+    :param pareto_front: List of individuals in the Pareto front
+    :return: The best individual based on TOPSIS ranking
+    """
+    # Step 1: Normalize the objective functions (cost and energy consumption)
+    min_lcc = min([ind.individual['lcc'] for ind in pareto_front])
+    max_lcc = max([ind.individual['lcc'] for ind in pareto_front])
+    min_lce = min([ind.individual['total_energy_consumption'] for ind in pareto_front])
+    max_lce = max([ind.individual['total_energy_consumption'] for ind in pareto_front])
+
+    normalized_pareto_front = []
+    for ind in pareto_front:
+      normalized_lcc = (ind.individual['lcc'] - min_lcc) / (max_lcc - min_lcc) if max_lcc > min_lcc else 0
+      normalized_lce = (ind.individual['total_energy_consumption'] - min_lce) / (
+            max_lce - min_lce) if max_lce > min_lce else 0
+      normalized_pareto_front.append((ind, normalized_lcc, normalized_lce))
+
+    # Step 2: Calculate the ideal and worst solutions
+    ideal_solution = (0, 0)  # Ideal is minimum LCC and minimum LCE (0, 0 after normalization)
+    worst_solution = (1, 1)  # Worst is maximum LCC and maximum LCE (1, 1 after normalization)
+
+    # Step 3: Calculate the distance to the ideal and worst solutions
+    best_distances = []
+    worst_distances = []
+
+    for ind, normalized_lcc, normalized_lce in normalized_pareto_front:
+      distance_to_ideal = math.sqrt(
+        (normalized_lcc - ideal_solution[0]) ** 2 + (normalized_lce - ideal_solution[1]) ** 2)
+      distance_to_worst = math.sqrt(
+        (normalized_lcc - worst_solution[0]) ** 2 + (normalized_lce - worst_solution[1]) ** 2)
+      best_distances.append(distance_to_ideal)
+      worst_distances.append(distance_to_worst)
+
+    # Step 4: Calculate relative closeness to the ideal solution
+    similarity = [worst / (best + worst) for best, worst in zip(best_distances, worst_distances)]
+
+    # Step 5: Select the individual with the highest similarity score
+    best_index = similarity.index(max(similarity))
+    best_solution = pareto_front[best_index]
+
+    return best_solution
+
+  @staticmethod
+  def design_period_identification(building):
+    def get_start_end_indices(max_day_index, total_days):
+      if max_day_index > 0 and max_day_index < total_days - 1:
+        start_index = (max_day_index - 1) * 24
+        end_index = (max_day_index + 2) * 24
+      elif max_day_index == 0:
+        start_index = 0
+        end_index = (max_day_index + 2) * 24
+      else:
+        start_index = (max_day_index - 1) * 24
+        end_index = total_days * 24
+      return start_index, end_index
+
+    # Calculate daily demands
+    heating_daily_demands = [sum(building.heating_demand[cte.HOUR][i:i + 24]) for i in
+                             range(0, len(building.heating_demand[cte.HOUR]), 24)]
+    cooling_daily_demands = [sum(building.cooling_demand[cte.HOUR][i:i + 24]) for i in
+                             range(0, len(building.cooling_demand[cte.HOUR]), 24)]
+    dhw_daily_demands = [sum(building.domestic_hot_water_heat_demand[cte.HOUR][i:i + 24]) for i in
+                         range(0, len(building.domestic_hot_water_heat_demand[cte.HOUR]), 24)]
+    # Get the day with maximum demand for each type
+    heating_max_day = heating_daily_demands.index(max(heating_daily_demands))
+    cooling_max_day = cooling_daily_demands.index(max(cooling_daily_demands))
+    dhw_max_day = dhw_daily_demands.index(max(dhw_daily_demands))
+    # Get the start and end indices for each demand type
+    heating_start, heating_end = get_start_end_indices(heating_max_day, len(heating_daily_demands))
+    cooling_start, cooling_end = get_start_end_indices(cooling_max_day, len(cooling_daily_demands))
+    dhw_start, dhw_end = get_start_end_indices(dhw_max_day, len(dhw_daily_demands))
+    # Return the design period energy demands
+    return {
+      f'{cte.HEATING}': {'demands': building.heating_demand[cte.HOUR][heating_start:heating_end],
+                         'start_index': heating_start, 'end_index': heating_end},
+      f'{cte.COOLING}': {'demands': building.cooling_demand[cte.HOUR][cooling_start:cooling_end],
+                         'start_index': cooling_start, 'end_index': cooling_end},
+      f'{cte.DOMESTIC_HOT_WATER}': {'demands': building.domestic_hot_water_heat_demand[cte.HOUR][dhw_start:dhw_end],
+                                    'start_index': dhw_start, 'end_index': dhw_end}
+    }
+

energy_system_modelling_package/energy_system_modelling_factories/system_sizing_methods/optimal_sizing.py

@@ -0,0 +1,42 @@
+import hub.helpers.constants as cte
+from energy_system_modelling_package.energy_system_modelling_factories.system_sizing_methods.genetic_algorithm.single_objective_genetic_algorithm import \
+  SingleObjectiveGeneticAlgorithm
+
+
+class OptimalSizing:
+  def __init__(self, city, optimization_scenario):
+    self.city = city
+    self.optimization_scenario = optimization_scenario
+
+  def enrich_buildings(self):
+    for building in self.city.buildings:
+      for energy_system in building.energy_systems:
+        if len(energy_system.generation_systems) == 1 and energy_system.generation_systems[0].energy_storage_systems is None:
+          if energy_system.generation_systems[0].system_type == cte.PHOTOVOLTAIC:
+            pass
+          else:
+            if cte.HEATING in energy_system.demand_types:
+              if cte.DOMESTIC_HOT_WATER in energy_system.demand_types:
+                design_load = max([building.heating_demand[cte.HOUR][i] +
+                                   building.domestic_hot_water_heat_demand[cte.HOUR][i] for i in
+                                   range(len(building.heating_demand))]) / cte.WATTS_HOUR_TO_JULES
+              else:
+                design_load = building.heating_peak_load[cte.YEAR][0]
+              energy_system.generation_systems[0].nominal_heat_output = design_load
+            elif cte.COOLING in energy_system.demand_types:
+              energy_system.generation_systems[0].nominal_cooling_output = building.cooling_peak_load[cte.YEAR][0]
+        else:
+          optimized_system = SingleObjectiveGeneticAlgorithm(optimization_scenario=self.optimization_scenario).solve_ga(building, energy_system)
+          for generation_system in energy_system.generation_systems:
+            system_type = generation_system.system_type
+            for generation_component in optimized_system.individual['Generation Components']:
+              if generation_component['type'] == system_type:
+                generation_system.nominal_heat_output = generation_component['heating_capacity']
+                generation_system.nominal_cooling_output = generation_component['cooling_capacity']
+            if generation_system.energy_storage_systems is not None:
+              for storage_system in generation_system.energy_storage_systems:
+                storage_type = f'{storage_system.type_energy_stored}_storage'
+                for storage_component in optimized_system.individual['Energy Storage Components']:
+                  if storage_component['type'] == storage_type:
+                    storage_system.nominal_capacity = storage_component['capacity']
+                    storage_system.volume = storage_component['volume']

energy_system_modelling_package/energy_system_modelling_factories/system_sizing_methods/peak_load_sizing.py

@@ -0,0 +1,99 @@
+import hub.helpers.constants as cte
+
+
+class PeakLoadSizing:
+  def __init__(self, city, default_primary_unit_percentage=0.7, storage_peak_coverage=3):
+    self.city = city
+    self.default_primary_unit_percentage = default_primary_unit_percentage
+    self.storage_peak_coverage = storage_peak_coverage
+
+  def enrich_buildings(self):
+    total_demand = 0
+    for building in self.city.buildings:
+      energy_systems = building.energy_systems
+      for energy_system in energy_systems:
+        if cte.HEATING in energy_system.demand_types:
+          if cte.DOMESTIC_HOT_WATER in energy_system.demand_types:
+            total_demand = [(building.heating_demand[cte.HOUR][i] +
+                             building.domestic_hot_water_heat_demand[cte.HOUR][i]) / cte.WATTS_HOUR_TO_JULES
+                            for i in range(len(building.heating_demand[cte.HOUR]))]
+          else:
+            total_demand = building.heating_peak_load[cte.YEAR]
+          design_load = max(total_demand)
+          self.allocate_capacity(energy_system, design_load, cte.HEATING, self.default_primary_unit_percentage)
+          if cte.COOLING in energy_system.demand_types:
+            cooling_design_load = building.cooling_peak_load[cte.YEAR][0]
+            self.allocate_capacity(energy_system, cooling_design_load, cte.COOLING,
+                                   self.default_primary_unit_percentage)
+
+        elif cte.COOLING in energy_system.demand_types:
+          design_load = building.cooling_peak_load[cte.YEAR][0]
+          self.allocate_capacity(energy_system, design_load, cte.COOLING, self.default_primary_unit_percentage)
+        elif cte.DOMESTIC_HOT_WATER in energy_system.demand_types:
+          design_load = building.domestic_hot_water_peak_load[cte.YEAR][0]
+          self.allocate_capacity(energy_system, design_load, cte.DOMESTIC_HOT_WATER,
+                                 self.default_primary_unit_percentage)
+
+        for generation_system in energy_system.generation_systems:
+          storage_systems = generation_system.energy_storage_systems
+          if storage_systems is not None:
+            if cte.DOMESTIC_HOT_WATER in energy_system.demand_types:
+              operation_range = 10
+            else:
+              operation_range = 20
+            for storage_system in storage_systems:
+              if storage_system.type_energy_stored == cte.THERMAL:
+                self.tes_sizing(storage_system, max(total_demand), self.storage_peak_coverage, operation_range)
+
+  def allocate_capacity(self, energy_system, design_load, demand_type, default_primary_unit_percentage):
+    if len(energy_system.generation_systems) == 1:
+      # If there's only one generation system, it gets the full design load.
+      if demand_type == cte.HEATING or demand_type == cte.DOMESTIC_HOT_WATER:
+        energy_system.generation_systems[0].nominal_heat_output = design_load
+      elif demand_type == cte.COOLING:
+        energy_system.generation_systems[0].nominal_cooling_output = design_load
+    else:
+      cooling_equipments_number = 0
+      # Distribute the load among generation systems.
+      max_efficiency = 0
+      main_generation_unit = None
+      for generation_system in energy_system.generation_systems:
+        if demand_type == cte.HEATING or demand_type == cte.DOMESTIC_HOT_WATER:
+          if max_efficiency < float(generation_system.heat_efficiency):
+            max_efficiency = float(generation_system.heat_efficiency)
+            main_generation_unit = generation_system
+        elif demand_type == cte.COOLING and generation_system.fuel_type == cte.ELECTRICITY:
+          cooling_equipments_number += 1
+          if max_efficiency < float(generation_system.cooling_efficiency):
+            max_efficiency = float(generation_system.heat_efficiency)
+            main_generation_unit = generation_system
+
+      for generation_system in energy_system.generation_systems:
+        if generation_system.system_type == main_generation_unit.system_type:
+          if demand_type == cte.HEATING or demand_type == cte.DOMESTIC_HOT_WATER:
+            generation_system.nominal_heat_output = round(default_primary_unit_percentage * design_load)
+          elif demand_type == cte.COOLING and cooling_equipments_number > 1:
+            generation_system.nominal_cooling_output = round(default_primary_unit_percentage * design_load)
+          else:
+            generation_system.nominal_cooling_output = design_load
+        else:
+          if demand_type == cte.HEATING or demand_type == cte.DOMESTIC_HOT_WATER:
+            generation_system.nominal_heat_output = round(((1 - default_primary_unit_percentage) * design_load /
+                                                     (len(energy_system.generation_systems) - 1)))
+          elif demand_type == cte.COOLING and cooling_equipments_number > 1:
+            generation_system.nominal_cooling_output = round(((1 - default_primary_unit_percentage) * design_load /
+                                                        (len(energy_system.generation_systems) - 1)))
+
+
+  @staticmethod
+  def tes_sizing(storage, peak_load, coverage, operational_temperature_range):
+    storage.volume = round((peak_load * coverage * cte.WATTS_HOUR_TO_JULES) /
+                           (cte.WATER_HEAT_CAPACITY * cte.WATER_DENSITY * operational_temperature_range))
+
+  @staticmethod
+  def cooling_equipments(energy_system):
+    counter = 0
+    for generation_system in energy_system.generation_systems:
+      if generation_system.fuel_type == cte.ELECTRICITY:
+        counter += 1
+    return counter

energy_system_modelling_package/energy_system_retrofit/energy_system_retrofit_report.py

@@ -0,0 +1,595 @@
+import os
+import hub.helpers.constants as cte
+import matplotlib.pyplot as plt
+from matplotlib import cm
+from energy_system_modelling_package.report_creation import LatexReport
+from matplotlib.ticker import MaxNLocator
+import numpy as np
+from pathlib import Path
+import glob
+
+
+class EnergySystemRetrofitReport:
+  def __init__(self, city, output_path, retrofit_scenario, current_status_energy_consumption_data,
+               retrofitted_energy_consumption_data, current_status_lcc_data, retrofitted_lcc_data):
+    self.city = city
+    self.current_status_data = current_status_energy_consumption_data
+    self.retrofitted_data = retrofitted_energy_consumption_data
+    self.current_status_lcc = current_status_lcc_data
+    self.retrofitted_lcc = retrofitted_lcc_data
+    self.output_path = output_path
+    self.content = []
+    self.retrofit_scenario = retrofit_scenario
+    self.report = LatexReport('energy_system_retrofit_report',
+                              'Energy System Retrofit Report', self.retrofit_scenario, output_path)
+    self.system_schemas_path = (Path(__file__).parent.parent.parent / 'hub' / 'data' / 'energy_systems' / 'schemas')
+    self.charts_path = Path(output_path) / 'charts'
+    self.charts_path.mkdir(parents=True, exist_ok=True)
+    files = glob.glob(f'{self.charts_path}/*')
+    for file in files:
+      os.remove(file)
+
+  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 $m^2$", "Heating Demand Intensity kWh/ $m^2$",
+                             "Cooling Demand Intensity kWh/ $m^2$", "Electricity Intensity kWh/ $m^2$"]]
+    peak_load_data = [["Building Name", "Heating Peak Load (kW)", "Cooling Peak Load (kW)",
+                       "Domestic Hot Water Peak Load (kW)"]]
+
+    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] / 3.6e9, '.2f')),
+        str(format(building.cooling_demand[cte.YEAR][0] / 3.6e9, '.2f')),
+        str(format(building.domestic_hot_water_heat_demand[cte.YEAR][0] / 3.6e9, '.2f')),
+        str(format(
+          (building.lighting_electrical_demand[cte.YEAR][0] + building.appliances_electrical_demand[cte.YEAR][0])
+          / 3.6e9, '.2f')),
+      ]
+      intensity_data = [
+        building.name,
+        str(format(total_floor_area, '.2f')),
+        str(format(building.heating_demand[cte.YEAR][0] / (3.6e6 * total_floor_area), '.2f')),
+        str(format(building.cooling_demand[cte.YEAR][0] / (3.6e6 * total_floor_area), '.2f')),
+        str(format(
+          (building.lighting_electrical_demand[cte.YEAR][0] + building.appliances_electrical_demand[cte.YEAR][0]) /
+          (3.6e6 * total_floor_area), '.2f'))
+      ]
+      peak_data = [
+        building.name,
+        str(format(building.heating_peak_load[cte.YEAR][0] / 1000, '.2f')),
+        str(format(building.cooling_peak_load[cte.YEAR][0] / 1000, '.2f')),
+        str(format(
+          (building.lighting_electrical_demand[cte.YEAR][0] + building.appliances_electrical_demand[cte.YEAR][0]) /
+          (3.6e6 * total_floor_area), '.2f'))
+      ]
+      table_data.append(building_data)
+      intensity_table_data.append(intensity_data)
+      peak_load_data.append(peak_data)
+
+    self.report.add_table(table_data, caption='Buildings Energy Consumption Data')
+    self.report.add_table(intensity_table_data, caption='Buildings Energy Use Intensity Data')
+    self.report.add_table(peak_load_data, caption='Buildings Peak Load Data')
+
+  def plot_monthly_energy_demands(self, data, file_name, title):
+    # Data preparation
+    months = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']
+    demands = {
+      'Heating': ('heating', '#2196f3'),
+      'Cooling': ('cooling', '#ff5a5f'),
+      'DHW': ('dhw', '#4caf50'),
+      'Electricity': ('lighting_appliance', '#ffc107')
+    }
+
+    # Helper function for plotting
+    def plot_bar_chart(ax, demand_type, color, ylabel, title):
+      values = data[demand_type]
+      ax.bar(months, values, color=color, width=0.6, zorder=2)
+      ax.grid(which="major", axis='x', color='#DAD8D7', alpha=0.5, zorder=1)
+      ax.grid(which="major", axis='y', color='#DAD8D7', alpha=0.5, zorder=1)
+      ax.set_ylabel(ylabel, fontsize=14, labelpad=10)
+      ax.set_title(title, fontsize=14, weight='bold', alpha=.8, pad=40)
+      ax.xaxis.set_major_locator(MaxNLocator(integer=True))
+      ax.yaxis.set_major_locator(MaxNLocator(integer=True))
+      ax.set_xticks(np.arange(len(months)))
+      ax.set_xticklabels(months, rotation=45, ha='right')
+      ax.bar_label(ax.containers[0], padding=3, color='black', fontsize=12, rotation=90)
+      ax.spines[['top', 'left', 'bottom']].set_visible(False)
+      ax.spines['right'].set_linewidth(1.1)
+      average_value = np.mean(values)
+      ax.axhline(y=average_value, color='grey', linewidth=2, linestyle='--')
+      ax.text(len(months) - 1, average_value, f'Average = {average_value:.1f} kWh', ha='right', va='bottom',
+              color='grey')
+
+    # Plotting
+    fig, axs = plt.subplots(4, 1, figsize=(20, 16), dpi=96)
+    fig.suptitle(title, fontsize=16, weight='bold', alpha=.8)
+
+    plot_bar_chart(axs[0], 'heating', demands['Heating'][1], 'Heating Demand (kWh)', 'Monthly Heating Demand')
+    plot_bar_chart(axs[1], 'cooling', demands['Cooling'][1], 'Cooling Demand (kWh)', 'Monthly Cooling Demand')
+    plot_bar_chart(axs[2], 'dhw', demands['DHW'][1], 'DHW Demand (kWh)', 'Monthly DHW Demand')
+    plot_bar_chart(axs[3], 'lighting_appliance', demands['Electricity'][1], 'Electricity Demand (kWh)',
+                   'Monthly Electricity Demand')
+
+    # Set a white background
+    fig.patch.set_facecolor('white')
+
+    # Adjust the margins around the plot area
+    plt.subplots_adjust(left=0.05, right=0.95, top=0.9, bottom=0.1, hspace=0.5)
+
+    # Save the plot
+    chart_path = self.charts_path / f'{file_name}.png'
+    plt.savefig(chart_path, bbox_inches='tight')
+    plt.close()
+
+    return chart_path
+
+  def plot_monthly_energy_consumption(self, data, file_name, title):
+    # Data preparation
+    months = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']
+    consumptions = {
+      'Heating': ('heating', '#2196f3', 'Heating Consumption (kWh)', 'Monthly Energy Consumption for Heating'),
+      'Cooling': ('cooling', '#ff5a5f', 'Cooling Consumption (kWh)', 'Monthly Energy Consumption for Cooling'),
+      'DHW': ('dhw', '#4caf50', 'DHW Consumption (kWh)', 'Monthly DHW Consumption')
+    }
+
+    # Helper function for plotting
+    def plot_bar_chart(ax, consumption_type, color, ylabel, title):
+      values = data[consumption_type]
+      ax.bar(months, values, color=color, width=0.6, zorder=2)
+      ax.grid(which="major", axis='x', color='#DAD8D7', alpha=0.5, zorder=1)
+      ax.grid(which="major", axis='y', color='#DAD8D7', alpha=0.5, zorder=1)
+      ax.set_xlabel('Month', fontsize=12, labelpad=10)
+      ax.set_ylabel(ylabel, fontsize=14, labelpad=10)
+      ax.set_title(title, fontsize=14, weight='bold', alpha=.8, pad=40)
+      ax.xaxis.set_major_locator(MaxNLocator(integer=True))
+      ax.yaxis.set_major_locator(MaxNLocator(integer=True))
+      ax.set_xticks(np.arange(len(months)))
+      ax.set_xticklabels(months, rotation=45, ha='right')
+      ax.bar_label(ax.containers[0], padding=3, color='black', fontsize=12, rotation=90)
+      ax.spines[['top', 'left', 'bottom']].set_visible(False)
+      ax.spines['right'].set_linewidth(1.1)
+      average_value = np.mean(values)
+      ax.axhline(y=average_value, color='grey', linewidth=2, linestyle='--')
+      ax.text(len(months) - 1, average_value, f'Average = {average_value:.1f} kWh', ha='right', va='bottom',
+              color='grey')
+
+    # Plotting
+    fig, axs = plt.subplots(3, 1, figsize=(20, 15), dpi=96)
+    fig.suptitle(title, fontsize=16, weight='bold', alpha=.8)
+
+    plot_bar_chart(axs[0], 'heating', consumptions['Heating'][1], consumptions['Heating'][2],
+                   consumptions['Heating'][3])
+    plot_bar_chart(axs[1], 'cooling', consumptions['Cooling'][1], consumptions['Cooling'][2],
+                   consumptions['Cooling'][3])
+    plot_bar_chart(axs[2], 'dhw', consumptions['DHW'][1], consumptions['DHW'][2], consumptions['DHW'][3])
+
+    # Set a white background
+    fig.patch.set_facecolor('white')
+
+    # Adjust the margins around the plot area
+    plt.subplots_adjust(left=0.05, right=0.95, top=0.9, bottom=0.1, wspace=0.3, hspace=0.5)
+
+    # Save the plot
+    chart_path = self.charts_path / f'{file_name}.png'
+    plt.savefig(chart_path, bbox_inches='tight')
+    plt.close()
+
+    return chart_path
+
+  def fuel_consumption_breakdown(self, file_name, data):
+    fuel_consumption_breakdown = {}
+    for building in self.city.buildings:
+      for key, breakdown in data[f'{building.name}']['energy_consumption_breakdown'].items():
+        if key not in fuel_consumption_breakdown:
+          fuel_consumption_breakdown[key] = {sector: 0 for sector in breakdown}
+        for sector, value in breakdown.items():
+          if sector in fuel_consumption_breakdown[key]:
+            fuel_consumption_breakdown[key][sector] += value / 3.6e6
+          else:
+            fuel_consumption_breakdown[key][sector] = value / 3.6e6
+
+    plt.style.use('ggplot')
+    num_keys = len(fuel_consumption_breakdown)
+    fig, axs = plt.subplots(1 if num_keys <= 2 else num_keys, min(num_keys, 2), figsize=(12, 5))
+    axs = axs if num_keys > 1 else [axs]  # Ensure axs is always iterable
+
+    for i, (fuel, breakdown) in enumerate(fuel_consumption_breakdown.items()):
+      labels = breakdown.keys()
+      values = breakdown.values()
+      colors = cm.get_cmap('tab20c', len(labels))
+      ax = axs[i] if num_keys > 1 else axs[0]
+      ax.pie(values, labels=labels,
+             autopct=lambda pct: f"{pct:.1f}%\n({pct / 100 * sum(values):.2f})",
+             startangle=90, colors=[colors(j) for j in range(len(labels))])
+      ax.set_title(f'{fuel} Consumption Breakdown')
+
+    plt.suptitle('City Energy Consumption Breakdown', fontsize=16, fontweight='bold')
+    plt.tight_layout(rect=[0, 0, 1, 0.95])  # Adjust layout to fit the suptitle
+
+    chart_path = self.charts_path / f'{file_name}.png'
+    plt.savefig(chart_path, dpi=300)
+    plt.close()
+    return chart_path
+
+  def energy_system_archetype_schematic(self):
+    energy_system_archetypes = {}
+    for building in self.city.buildings:
+      if building.energy_systems_archetype_name not in energy_system_archetypes:
+        energy_system_archetypes[building.energy_systems_archetype_name] = [building.name]
+      else:
+        energy_system_archetypes[building.energy_systems_archetype_name].append(building.name)
+
+    text = ""
+    items = ""
+    for archetype, buildings in energy_system_archetypes.items():
+      buildings_str = ", ".join(buildings)
+      text += f"Figure 4 shows the schematic of the proposed energy system for buildings {buildings_str}.\n"
+      if archetype in ['PV+4Pipe+DHW', 'PV+ASHP+GasBoiler+TES', 'Central 4 Pipes Air to Water Heat Pump and Gas Boiler with Independent Water Heating and PV']:
+        text += "This energy system archetype is formed of the following systems: \par"
+        items = ['Rooftop Photovoltaic System: The rooftop PV system is tied to the grid and in case there is surplus '
+                 'energy, it is sold to Hydro-Quebec through their Net-Meterin program.',
+                 '4-Pipe HVAC System: This systems includes a 4-pipe heat pump capable of generating heat and cooling '
+                 'at the same time, a natural gas boiler as the auxiliary heating system, and a hot water storage tank.'
+                 'The temperature inside the tank is kept between 40-55 C. The cooling demand is totally supplied by '
+                 'the heat pump unit.',
+                 'Domestic Hot Water Heat Pump System: This system is in charge of supplying domestic hot water demand.'
+                 'The heat pump is connected to a thermal storage tank with electric resistance heating coil inside it.'
+                 ' The temperature inside the tank should always remain above 60 C.']
+
+      self.report.add_text(text)
+      self.report.add_itemize(items=items)
+      schema_path = self.system_schemas_path / f'{archetype}.jpg'
+      self.report.add_image(str(schema_path).replace('\\', '/'),
+                            f'Proposed energy system for buildings {buildings_str}')
+
+  def plot_monthly_radiation(self):
+    months = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']
+    monthly_roof_radiation = []
+    for i in range(len(months)):
+      tilted_radiation = 0
+      for building in self.city.buildings:
+        tilted_radiation += (building.roofs[0].global_irradiance_tilted[cte.MONTH][i] / 1000)
+      monthly_roof_radiation.append(tilted_radiation)
+
+    def plot_bar_chart(ax, months, values, color, ylabel, title):
+      ax.bar(months, values, color=color, width=0.6, zorder=2)
+      ax.grid(which="major", axis='x', color='#DAD8D7', alpha=0.5, zorder=1)
+      ax.grid(which="major", axis='y', color='#DAD8D7', alpha=0.5, zorder=1)
+      ax.set_xlabel('Month', fontsize=12, labelpad=10)
+      ax.set_ylabel(ylabel, fontsize=14, labelpad=10)
+      ax.set_title(title, fontsize=14, weight='bold', alpha=.8, pad=40)
+      ax.xaxis.set_major_locator(MaxNLocator(integer=True))
+      ax.yaxis.set_major_locator(MaxNLocator(integer=True))
+      ax.set_xticks(np.arange(len(months)))
+      ax.set_xticklabels(months, rotation=45, ha='right')
+      ax.bar_label(ax.containers[0], padding=3, color='black', fontsize=12, rotation=90)
+      ax.spines[['top', 'left', 'bottom']].set_visible(False)
+      ax.spines['right'].set_linewidth(1.1)
+      average_value = np.mean(values)
+      ax.axhline(y=average_value, color='grey', linewidth=2, linestyle='--')
+      ax.text(len(months) - 1, average_value, f'Average = {average_value:.1f} kWh', ha='right', va='bottom',
+              color='grey')
+
+    # Plotting the bar chart
+    fig, ax = plt.subplots(figsize=(15, 8), dpi=96)
+    plot_bar_chart(ax, months, monthly_roof_radiation, '#ffc107', 'Tilted Roof Radiation (kWh / m2)',
+                   'Monthly Tilted Roof Radiation')
+
+    # Set a white background
+    fig.patch.set_facecolor('white')
+
+    # Adjust the margins around the plot area
+    plt.subplots_adjust(left=0.1, right=0.95, top=0.9, bottom=0.1)
+
+    # Save the plot
+    chart_path = self.charts_path / 'monthly_tilted_roof_radiation.png'
+    plt.savefig(chart_path, bbox_inches='tight')
+    plt.close()
+    return chart_path
+
+  def energy_consumption_comparison(self, current_status_data, retrofitted_data, file_name, title):
+    months = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']
+    consumptions = {
+      'Heating': ('heating', '#2196f3', 'Heating Consumption (kWh)', 'Monthly Energy Consumption for Heating'),
+      'Cooling': ('cooling', '#ff5a5f', 'Cooling Consumption (kWh)', 'Monthly Energy Consumption for Cooling'),
+      'DHW': ('dhw', '#4caf50', 'DHW Consumption (kWh)', 'Monthly DHW Consumption')
+    }
+
+    # Helper function for plotting
+    def plot_double_bar_chart(ax, consumption_type, color, ylabel, title):
+      current_values = current_status_data[consumption_type]
+      retrofitted_values = retrofitted_data[consumption_type]
+      bar_width = 0.35
+      index = np.arange(len(months))
+
+      ax.bar(index, current_values, bar_width, label='Current Status', color=color, alpha=0.7, zorder=2)
+      ax.bar(index + bar_width, retrofitted_values, bar_width, label='Retrofitted', color=color, hatch='/', zorder=2)
+
+      ax.grid(which="major", axis='x', color='#DAD8D7', alpha=0.5, zorder=1)
+      ax.grid(which="major", axis='y', color='#DAD8D7', alpha=0.5, zorder=1)
+      ax.set_xlabel('Month', fontsize=12, labelpad=10)
+      ax.set_ylabel(ylabel, fontsize=14, labelpad=10)
+      ax.set_title(title, fontsize=14, weight='bold', alpha=.8, pad=40)
+      ax.set_xticks(index + bar_width / 2)
+      ax.set_xticklabels(months, rotation=45, ha='right')
+      ax.legend()
+
+      # Adding bar labels
+      ax.bar_label(ax.containers[0], padding=3, color='black', fontsize=12, rotation=90)
+      ax.bar_label(ax.containers[1], padding=3, color='black', fontsize=12, rotation=90)
+
+      ax.spines[['top', 'left', 'bottom']].set_visible(False)
+      ax.spines['right'].set_linewidth(1.1)
+
+    # Plotting
+    fig, axs = plt.subplots(3, 1, figsize=(20, 25), dpi=96)
+    fig.suptitle(title, fontsize=16, weight='bold', alpha=.8)
+
+    plot_double_bar_chart(axs[0], 'heating', consumptions['Heating'][1], consumptions['Heating'][2],
+                          consumptions['Heating'][3])
+    plot_double_bar_chart(axs[1], 'cooling', consumptions['Cooling'][1], consumptions['Cooling'][2],
+                          consumptions['Cooling'][3])
+    plot_double_bar_chart(axs[2], 'dhw', consumptions['DHW'][1], consumptions['DHW'][2], consumptions['DHW'][3])
+
+    # Set a white background
+    fig.patch.set_facecolor('white')
+
+    # Adjust the margins around the plot area
+    plt.subplots_adjust(left=0.05, right=0.95, top=0.9, bottom=0.1, wspace=0.3, hspace=0.5)
+
+    # Save the plot
+    chart_path = self.charts_path / f'{file_name}.png'
+    plt.savefig(chart_path, bbox_inches='tight')
+    plt.close()
+
+    return chart_path
+
+  def yearly_consumption_comparison(self):
+    current_total_consumption = round(self.current_status_data['total_consumption'], 2)
+    retrofitted_total_consumption = round(self.retrofitted_data['total_consumption'], 2)
+    text = (
+      f'The total yearly energy consumption before and after the retrofit are {current_total_consumption} MWh and '
+      f'{retrofitted_total_consumption} MWh, respectively.')
+    if retrofitted_total_consumption < current_total_consumption:
+      change = str(round((current_total_consumption - retrofitted_total_consumption) * 100 / current_total_consumption,
+                         2))
+      text += f'Therefore, the total yearly energy consumption decreased by {change} \%.'
+    else:
+      change = str(round((retrofitted_total_consumption - current_total_consumption) * 100 /
+                         retrofitted_total_consumption, 2))
+      text += f'Therefore, the total yearly energy consumption increased by {change} \%. \par'
+    self.report.add_text(text)
+
+  def pv_system(self):
+    self.report.add_text('The first step in PV assessments is evaluating the potential of buildings for installing '
+                         'rooftop PV system. The benchmark value used for this evaluation is the mean yearly solar '
+                         'incident in Montreal. According to Hydro-Quebec, the mean annual incident in Montreal is 1350'
+                         'kWh/m2. Therefore, any building with rooftop annual global horizontal radiation of less than '
+                         '1080 kWh/m2 is considered to be infeasible. Table 4 shows the yearly horizontal radiation on '
+                         'buildings roofs. \par')
+    radiation_data = [
+      ["Building Name", "Roof Area $m^2$", "Function", "Rooftop Annual Global Horizontal Radiation kWh/ $m^2$"]
+    ]
+    pv_feasible_buildings = []
+    for building in self.city.buildings:
+      if building.roofs[0].global_irradiance[cte.YEAR][0] > 1080:
+        pv_feasible_buildings.append(building.name)
+      data = [building.name, str(format(building.roofs[0].perimeter_area, '.2f')), building.function,
+              str(format(building.roofs[0].global_irradiance[cte.YEAR][0] / (cte.WATTS_HOUR_TO_JULES * 1000), '.2f'))]
+      radiation_data.append(data)
+
+    self.report.add_table(radiation_data,
+                          caption='Buildings Roof Characteristics')
+
+    if len(pv_feasible_buildings) == len(self.city.buildings):
+      buildings_str = 'all'
+    else:
+      buildings_str = ", ".join(pv_feasible_buildings)
+    self.report.add_text(f"From the table it can be seen that {buildings_str} buildings are good candidates to have "
+                         f"rooftop PV system. The next step is calculating the amount of solar radiation on a tilted "
+                         f"surface. Figure 5 shows the total monthly solar radiation on a surface with the tilt angle "
+                         f"of 45 degrees on the roofs of those buildings that are identified to have rooftop PV system."
+                         f"\par")
+    tilted_radiation = self.plot_monthly_radiation()
+    self.report.add_image(str(tilted_radiation).replace('\\', '/'),
+                          caption='Total Monthly Solar Radiation on Buildings Roofs on a 45 Degrees Tilted Surface',
+                          placement='H')
+    self.report.add_text('The first step in sizing the PV system is to find the available roof area. '
+                         'Few considerations need to be made here. The considerations include space for maintenance '
+                         'crew, space for mechanical equipment, and orientation correction factor to make sure all '
+                         'the panel are truly facing south. After all these considerations, the minimum distance '
+                         'between the panels to avoid shading throughout the year is found. Table 5 shows the number of'
+                         'panles on each buildings roof, yearly PV production, total electricity consumption, and self '
+                         'consumption. \par')
+
+    pv_output_table = [['Building Name', 'Total Surface Area of PV Panels ($m^2$)',
+                        'Total Solar Incident on PV Modules (MWh)', 'Yearly PV production (MWh)']]
+
+    for building in self.city.buildings:
+      if building.name in pv_feasible_buildings:
+        pv_data = []
+        pv_data.append(building.name)
+        yearly_solar_incident = (building.roofs[0].global_irradiance_tilted[cte.YEAR][0] *
+                                 building.roofs[0].installed_solar_collector_area) / (cte.WATTS_HOUR_TO_JULES * 1e6)
+        yearly_solar_incident_str = format(yearly_solar_incident, '.2f')
+        yearly_pv_output = building.onsite_electrical_production[cte.YEAR][0] / (cte.WATTS_HOUR_TO_JULES * 1e6)
+        yearly_pv_output_str = format(yearly_pv_output, '.2f')
+
+        pv_data.append(format(building.roofs[0].installed_solar_collector_area, '.2f'))
+        pv_data.append(yearly_solar_incident_str)
+        pv_data.append(yearly_pv_output_str)
+
+        pv_output_table.append(pv_data)
+
+    self.report.add_table(pv_output_table, caption='PV System Simulation Results', first_column_width=3)
+
+  def life_cycle_cost_stacked_bar(self, file_name, title):
+    # Aggregate LCC components for current and retrofitted statuses
+    current_status_capex = 0
+    current_status_opex = 0
+    current_status_maintenance = 0
+    current_status_end_of_life = 0
+    retrofitted_capex = 0
+    retrofitted_opex = 0
+    retrofitted_maintenance = 0
+    retrofitted_end_of_life = 0
+    current_status_operational_income = 0
+    retrofitted_operational_income = 0
+
+    for building in self.city.buildings:
+      current_status_capex += self.current_status_lcc[f'{building.name}']['capital_cost_per_sqm']
+      retrofitted_capex += self.retrofitted_lcc[f'{building.name}']['capital_cost_per_sqm']
+      current_status_opex += self.current_status_lcc[f'{building.name}']['operational_cost_per_sqm']
+      retrofitted_opex += self.retrofitted_lcc[f'{building.name}']['operational_cost_per_sqm']
+      current_status_maintenance += self.current_status_lcc[f'{building.name}']['maintenance_cost_per_sqm']
+      retrofitted_maintenance += self.retrofitted_lcc[f'{building.name}']['maintenance_cost_per_sqm']
+      current_status_end_of_life += self.current_status_lcc[f'{building.name}']['end_of_life_cost_per_sqm']
+      retrofitted_end_of_life += self.retrofitted_lcc[f'{building.name}']['end_of_life_cost_per_sqm']
+      current_status_operational_income += self.current_status_lcc[f'{building.name}']['operational_income_per_sqm']
+      retrofitted_operational_income += self.retrofitted_lcc[f'{building.name}']['operational_income_per_sqm']
+
+    current_status_lcc_components_sqm = {
+      'Capital Cost': current_status_capex / len(self.city.buildings),
+      'Operational Cost': (current_status_opex - current_status_operational_income) / len(self.city.buildings),
+      'Maintenance Cost': current_status_maintenance / len(self.city.buildings),
+      'End of Life Cost': current_status_end_of_life / len(self.city.buildings),
+    }
+    retrofitted_lcc_components_sqm = {
+      'Capital Cost': retrofitted_capex / len(self.city.buildings),
+      'Operational Cost': (retrofitted_opex - retrofitted_operational_income) / len(self.city.buildings),
+      'Maintenance Cost': retrofitted_maintenance / len(self.city.buildings),
+      'End of Life Cost': retrofitted_end_of_life / len(self.city.buildings),
+    }
+
+    labels = ['Current Status', 'Retrofitted Status']
+    categories = ['Capital Cost', 'Operational Cost', 'Maintenance Cost', 'End of Life Cost']
+    colors = ['#2196f3', '#ff5a5f', '#4caf50', '#ffc107']  # Added new color
+
+    # Data preparation
+    bar_width = 0.35
+    r = np.arange(len(labels))
+
+    fig, ax = plt.subplots(figsize=(12, 8), dpi=96)
+    fig.suptitle(title, fontsize=16, weight='bold', alpha=.8)
+
+    # Plotting current status data
+    bottom = np.zeros(len(labels))
+    for category, color in zip(categories, colors):
+      values = [current_status_lcc_components_sqm[category], retrofitted_lcc_components_sqm[category]]
+      ax.bar(r, values, bottom=bottom, color=color, edgecolor='white', width=bar_width, label=category)
+      bottom += values
+
+    # Adding summation annotations at the top of the bars
+    for idx, (x, total) in enumerate(zip(r, bottom)):
+      ax.text(x, total, f'{total:.1f}', ha='center', va='bottom', fontsize=12, fontweight='bold')
+
+    # Adding labels, title, and grid
+    ax.set_xlabel('LCC Components', fontsize=12, labelpad=10)
+    ax.set_ylabel('Average Cost (CAD/m²)', fontsize=14, labelpad=10)
+    ax.grid(which="major", axis='y', color='#DAD8D7', alpha=0.5, zorder=1)
+    ax.set_xticks(r)
+    ax.set_xticklabels(labels, rotation=45, ha='right')
+    ax.legend()
+
+    # Adding a white background
+    fig.patch.set_facecolor('white')
+
+    # Adjusting the margins around the plot area
+    plt.subplots_adjust(left=0.05, right=0.95, top=0.9, bottom=0.2)
+
+    # Save the plot
+    chart_path = self.charts_path / f'{file_name}.png'
+    plt.savefig(chart_path, bbox_inches='tight')
+    plt.close()
+
+    return chart_path
+
+  def create_report(self):
+    # Add sections and text to the report
+    self.report.add_section('Overview of the Current Status in Buildings')
+    self.report.add_text('In this section, an overview of the current status of buildings characteristics, '
+                         'energy demand and consumptions are provided')
+    self.report.add_subsection('Buildings Characteristics')
+
+    self.building_energy_info()
+
+    # current monthly demands and consumptions
+    current_monthly_demands = self.current_status_data['monthly_demands']
+    current_monthly_consumptions = self.current_status_data['monthly_consumptions']
+
+    # Plot and save demand chart
+    current_demand_chart_path = self.plot_monthly_energy_demands(data=current_monthly_demands,
+                                                                 file_name='current_monthly_demands',
+                                                                 title='Current Status Monthly Energy Demands')
+    # Plot and save consumption chart
+    current_consumption_chart_path = self.plot_monthly_energy_consumption(data=current_monthly_consumptions,
+                                                                          file_name='monthly_consumptions',
+                                                                          title='Monthly Energy Consumptions')
+    current_consumption_breakdown_path = self.fuel_consumption_breakdown('City_Energy_Consumption_Breakdown',
+                                                                         self.current_status_data)
+    retrofitted_consumption_breakdown_path = self.fuel_consumption_breakdown(
+      'fuel_consumption_breakdown_after_retrofit',
+      self.retrofitted_data)
+    life_cycle_cost_sqm_stacked_bar_chart_path = self.life_cycle_cost_stacked_bar('lcc_per_sqm',
+                                                                                  'LCC Analysis')
+    # Add current state of energy demands in the city
+    self.report.add_subsection('Current State of Energy Demands in the City')
+    self.report.add_text('The total monthly energy demands in the city are shown in Figure 1. It should be noted '
+                         'that the electricity demand refers to total lighting and appliance electricity demands')
+    self.report.add_image(str(current_demand_chart_path).replace('\\', '/'),
+                          'Total Monthly Energy Demands in City',
+                          placement='h')
+
+    # Add current state of energy consumption in the city
+    self.report.add_subsection('Current State of Energy Consumption in the City')
+    self.report.add_text('The following figure shows the amount of energy consumed to supply heating, cooling, and '
+                         'domestic hot water needs in the city. The details of the systems in each building before '
+                         'and after retrofit are provided in Section 4. \par')
+    self.report.add_image(str(current_consumption_chart_path).replace('\\', '/'),
+                          'Total Monthly Energy Consumptions in City',
+                          placement='H')
+    self.report.add_text('Figure 3 shows the yearly energy supplied to the city by fuel in different sectors. '
+                         'All the values are in kWh.')
+    self.report.add_image(str(current_consumption_breakdown_path).replace('\\', '/'),
+                          'Current Energy Consumption Breakdown in the City by Fuel',
+                          placement='H')
+    self.report.add_section(f'{self.retrofit_scenario}')
+    self.report.add_subsection('Proposed Systems')
+    self.energy_system_archetype_schematic()
+    if 'PV' in self.retrofit_scenario:
+      self.report.add_subsection('Rooftop Photovoltaic System Implementation')
+      self.pv_system()
+    if 'System' in self.retrofit_scenario:
+      self.report.add_subsection('Retrofitted HVAC and DHW Systems')
+      self.report.add_text('Figure 6 shows a comparison between total monthly energy consumption in the selected '
+                           'buildings before and after retrofitting.')
+      consumption_comparison = self.energy_consumption_comparison(self.current_status_data['monthly_consumptions'],
+                                                                  self.retrofitted_data['monthly_consumptions'],
+                                                                  'energy_consumption_comparison_in_city',
+                                                                  'Total Monthly Energy Consumption Comparison in '
+                                                                  'Buildings')
+      self.report.add_image(str(consumption_comparison).replace('\\', '/'),
+                            caption='Comparison of Total Monthly Energy Consumption in City Buildings',
+                            placement='H')
+      self.yearly_consumption_comparison()
+      self.report.add_text('Figure 7 shows the fuel consumption breakdown in the area after the retrofit.')
+      self.report.add_image(str(retrofitted_consumption_breakdown_path).replace('\\', '/'),
+                            caption=f'Fuel Consumption Breakdown After {self.retrofit_scenario}',
+                            placement='H')
+    self.report.add_subsection('Life Cycle Cost Analysis')
+    self.report.add_image(str(life_cycle_cost_sqm_stacked_bar_chart_path).replace('\\', '/'),
+                          caption='Average Life Cycle Cost Components',
+                          placement='H')
+
+    # Save and compile the report
+    self.report.save_report()
+    self.report.compile_to_pdf()

energy_system_modelling_package/energy_system_retrofit/energy_system_retrofit_results.py

@@ -0,0 +1,176 @@
+import hub.helpers.constants as cte
+
+
+def hourly_electricity_consumption_profile(building):
+  hourly_electricity_consumption = []
+  energy_systems = building.energy_systems
+  appliance = building.appliances_electrical_demand[cte.HOUR]
+  lighting = building.lighting_electrical_demand[cte.HOUR]
+  elec_heating = 0
+  elec_cooling = 0
+  elec_dhw = 0
+  if cte.HEATING in building.energy_consumption_breakdown[cte.ELECTRICITY]:
+    elec_heating = 1
+  if cte.COOLING in building.energy_consumption_breakdown[cte.ELECTRICITY]:
+    elec_cooling = 1
+  if cte.DOMESTIC_HOT_WATER in building.energy_consumption_breakdown[cte.ELECTRICITY]:
+    elec_dhw = 1
+  heating = None
+  cooling = None
+  dhw = None
+  if elec_heating == 1:
+    for energy_system in energy_systems:
+      if cte.HEATING in energy_system.demand_types:
+        for generation_system in energy_system.generation_systems:
+          if generation_system.fuel_type == cte.ELECTRICITY:
+            if cte.HEATING in generation_system.energy_consumption:
+              heating = generation_system.energy_consumption[cte.HEATING][cte.HOUR]
+            else:
+              if len(energy_system.generation_systems) > 1:
+                heating = [x / 2 for x in building.heating_consumption[cte.HOUR]]
+              else:
+                heating = building.heating_consumption[cte.HOUR]
+  if elec_dhw == 1:
+    for energy_system in energy_systems:
+      if cte.DOMESTIC_HOT_WATER in energy_system.demand_types:
+        for generation_system in energy_system.generation_systems:
+          if generation_system.fuel_type == cte.ELECTRICITY:
+            if cte.DOMESTIC_HOT_WATER in generation_system.energy_consumption:
+              dhw = generation_system.energy_consumption[cte.DOMESTIC_HOT_WATER][cte.HOUR]
+            else:
+              if len(energy_system.generation_systems) > 1:
+                dhw = [x / 2 for x in building.domestic_hot_water_consumption[cte.HOUR]]
+              else:
+                dhw = building.domestic_hot_water_consumption[cte.HOUR]
+
+  if elec_cooling == 1:
+    for energy_system in energy_systems:
+      if cte.COOLING in energy_system.demand_types:
+        for generation_system in energy_system.generation_systems:
+          if cte.COOLING in generation_system.energy_consumption:
+            cooling = generation_system.energy_consumption[cte.COOLING][cte.HOUR]
+          else:
+            if len(energy_system.generation_systems) > 1:
+              cooling = [x / 2 for x in building.cooling_consumption[cte.HOUR]]
+            else:
+              cooling = building.cooling_consumption[cte.HOUR]
+
+  for i in range(len(building.heating_demand[cte.HOUR])):
+    hourly = 0
+    hourly += appliance[i] / 3600
+    hourly += lighting[i] / 3600
+    if heating is not None:
+      hourly += heating[i] / 3600
+    if cooling is not None:
+      hourly += cooling[i] / 3600
+    if dhw is not None:
+      hourly += dhw[i] / 3600
+    hourly_electricity_consumption.append(hourly)
+  return hourly_electricity_consumption
+
+
+def consumption_data(city):
+  energy_consumption_data = {}
+  for building in city.buildings:
+    hourly_electricity_consumption = hourly_electricity_consumption_profile(building)
+    energy_consumption_data[f'{building.name}'] = {'heating_consumption': building.heating_consumption,
+                                                   'cooling_consumption': building.cooling_consumption,
+                                                   'domestic_hot_water_consumption':
+                                                   building.domestic_hot_water_consumption,
+                                                   'energy_consumption_breakdown':
+                                                   building.energy_consumption_breakdown,
+                                                   'hourly_electricity_consumption': hourly_electricity_consumption}
+  peak_electricity_consumption = 0
+  for building in energy_consumption_data:
+    peak_electricity_consumption += max(energy_consumption_data[building]['hourly_electricity_consumption'])
+  heating_demand_monthly = []
+  cooling_demand_monthly = []
+  dhw_demand_monthly = []
+  lighting_appliance_monthly = []
+  heating_consumption_monthly = []
+  cooling_consumption_monthly = []
+  dhw_consumption_monthly = []
+  for i in range(12):
+    heating_demand = 0
+    cooling_demand = 0
+    dhw_demand = 0
+    lighting_appliance_demand = 0
+    heating_consumption = 0
+    cooling_consumption = 0
+    dhw_consumption = 0
+    for building in city.buildings:
+      heating_demand += building.heating_demand[cte.MONTH][i] / 3.6e6
+      cooling_demand += building.cooling_demand[cte.MONTH][i] / 3.6e6
+      dhw_demand += building.domestic_hot_water_heat_demand[cte.MONTH][i] / 3.6e6
+      lighting_appliance_demand += building.lighting_electrical_demand[cte.MONTH][i] / 3.6e6
+      heating_consumption += building.heating_consumption[cte.MONTH][i] / 3.6e6
+      if building.cooling_consumption[cte.YEAR][0] == 0:
+        cooling_consumption += building.cooling_demand[cte.MONTH][i] / (3.6e6 * 2)
+      else:
+        cooling_consumption += building.cooling_consumption[cte.MONTH][i] / 3.6e6
+      dhw_consumption += building.domestic_hot_water_consumption[cte.MONTH][i] / 3.6e6
+    heating_demand_monthly.append(heating_demand)
+    cooling_demand_monthly.append(cooling_demand)
+    dhw_demand_monthly.append(dhw_demand)
+    lighting_appliance_monthly.append(lighting_appliance_demand)
+    heating_consumption_monthly.append(heating_consumption)
+    cooling_consumption_monthly.append(cooling_consumption)
+    dhw_consumption_monthly.append(dhw_consumption)
+
+  monthly_demands = {'heating': heating_demand_monthly,
+                     'cooling': cooling_demand_monthly,
+                     'dhw': dhw_demand_monthly,
+                     'lighting_appliance': lighting_appliance_monthly}
+  monthly_consumptions = {'heating': heating_consumption_monthly,
+                          'cooling': cooling_consumption_monthly,
+                          'dhw': dhw_consumption_monthly}
+  yearly_heating = 0
+  yearly_cooling = 0
+  yearly_dhw = 0
+  yearly_appliance = 0
+  yearly_lighting = 0
+  for building in city.buildings:
+    yearly_appliance += building.appliances_electrical_demand[cte.YEAR][0] / 3.6e9
+    yearly_lighting += building.lighting_electrical_demand[cte.YEAR][0] / 3.6e9
+    yearly_heating += building.heating_consumption[cte.YEAR][0] / 3.6e9
+    yearly_cooling += building.cooling_consumption[cte.YEAR][0] / 3.6e9
+    yearly_dhw += building.domestic_hot_water_consumption[cte.YEAR][0] / 3.6e9
+
+  total_consumption = yearly_heating + yearly_cooling + yearly_dhw + yearly_appliance + yearly_lighting
+  energy_consumption_data['monthly_demands'] = monthly_demands
+  energy_consumption_data['monthly_consumptions'] = monthly_consumptions
+  energy_consumption_data['total_consumption'] = total_consumption
+  energy_consumption_data['maximum_hourly_electricity_consumption'] = peak_electricity_consumption
+
+  return energy_consumption_data
+
+
+def cost_data(building, lcc_dataframe, cost_retrofit_scenario):
+  total_floor_area = 0
+  for thermal_zone in building.thermal_zones_from_internal_zones:
+    total_floor_area += thermal_zone.total_floor_area
+  capital_cost = lcc_dataframe.loc['total_capital_costs_systems', f'Scenario {cost_retrofit_scenario}']
+  operational_cost = lcc_dataframe.loc['total_operational_costs', f'Scenario {cost_retrofit_scenario}']
+  maintenance_cost = lcc_dataframe.loc['total_maintenance_costs', f'Scenario {cost_retrofit_scenario}']
+  end_of_life_cost = lcc_dataframe.loc['end_of_life_costs', f'Scenario {cost_retrofit_scenario}']
+  operational_income = lcc_dataframe.loc['operational_incomes', f'Scenario {cost_retrofit_scenario}']
+  total_life_cycle_cost = capital_cost + operational_cost + maintenance_cost + end_of_life_cost + operational_income
+  specific_capital_cost = capital_cost / total_floor_area
+  specific_operational_cost = operational_cost / total_floor_area
+  specific_maintenance_cost = maintenance_cost / total_floor_area
+  specific_end_of_life_cost = end_of_life_cost / total_floor_area
+  specific_operational_income = operational_income / total_floor_area
+  specific_life_cycle_cost = total_life_cycle_cost / total_floor_area
+  life_cycle_cost_analysis = {'capital_cost': capital_cost,
+                              'capital_cost_per_sqm': specific_capital_cost,
+                              'operational_cost': operational_cost,
+                              'operational_cost_per_sqm': specific_operational_cost,
+                              'maintenance_cost': maintenance_cost,
+                              'maintenance_cost_per_sqm': specific_maintenance_cost,
+                              'end_of_life_cost': end_of_life_cost,
+                              'end_of_life_cost_per_sqm': specific_end_of_life_cost,
+                              'operational_income': operational_income,
+                              'operational_income_per_sqm': specific_operational_income,
+                              'total_life_cycle_cost': total_life_cycle_cost,
+                              'total_life_cycle_cost_per_sqm': specific_life_cycle_cost}
+  return life_cycle_cost_analysis

energy_system_modelling_package/random_assignation.py

@@ -0,0 +1,123 @@
+"""
+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': 15,
+                                  'system 1 electricity': 35,
+                                  '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': 15,
+                                  'system 8 electricity': 35}
+
+residential_new_systems_percentage = {
+    'Central Hydronic Air and Gas Source Heating System with Unitary Split Cooling and Air Source HP DHW and PV': 100,
+    'Central Hydronic Air and Electricity Source Heating System with Unitary Split Cooling and Air Source HP DHW and PV': 0,
+    'Central Hydronic Ground and Gas Source Heating System with Unitary Split Cooling and Air Source HP DHW and PV': 0,
+    'Central Hydronic Ground and Electricity Source Heating System with Unitary Split Cooling and Air Source HP DHW '
+    'and PV': 0,
+    'Central Hydronic Water and Gas Source Heating System with Unitary Split Cooling and Air Source HP DHW and PV': 0,
+    'Central Hydronic Water and Electricity Source Heating System with Unitary Split Cooling and Air Source HP DHW '
+    'and PV': 0,
+    'Central Hydronic Air and Gas Source Heating System with Unitary Split and Air Source HP DHW': 0,
+    'Central Hydronic Air and Electricity Source Heating System with Unitary Split and Air Source HP DHW': 0,
+    'Central Hydronic Ground and Gas Source Heating System with Unitary Split and Air Source HP DHW': 0,
+    'Central Hydronic Ground and Electricity Source Heating System with Unitary Split and Air Source HP DHW': 0,
+    'Central Hydronic Water and Gas Source Heating System with Unitary Split and Air Source HP DHW': 0,
+    'Central Hydronic Water and Electricity Source Heating System with Unitary Split and Air Source HP DHW': 0,
+    'Rooftop PV System': 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
+

energy_system_modelling_package/report_creation.py

@@ -0,0 +1,119 @@
+import subprocess
+import datetime
+import os
+from pathlib import Path
+
+
+class LatexReport:
+  def __init__(self, file_name, title, subtitle, output_path):
+    self.file_name = file_name
+    self.output_path = Path(output_path) / 'report'
+    self.output_path.mkdir(parents=True, exist_ok=True)
+    self.file_path = self.output_path / f"{file_name}.tex"
+    self.content = []
+    self.content.append(r'\documentclass{article}')
+    self.content.append(r'\usepackage[margin=2.5cm]{geometry}')
+    self.content.append(r'\usepackage{graphicx}')
+    self.content.append(r'\usepackage{tabularx}')
+    self.content.append(r'\usepackage{multirow}')
+    self.content.append(r'\usepackage{float}')
+    self.content.append(r'\begin{document}')
+
+    current_datetime = datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S")
+
+    self.content.append(r'\title{' + title + '}')
+    self.content.append(r'\author{Next-Generation Cities Institute}')
+    self.content.append(r'\date{}')
+    self.content.append(r'\maketitle')
+
+    self.content.append(r'\begin{center}')
+    self.content.append(r'\large ' + subtitle + r'\\')
+    self.content.append(r'\large ' + current_datetime)
+    self.content.append(r'\end{center}')
+
+  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_subsubsection(self, subsection_title):
+    self.content.append(r'\subsubsection{' + subsection_title + r'}')
+
+  def add_text(self, text):
+    self.content.append(text)
+
+  def add_table(self, table_data, caption=None, first_column_width=None, merge_first_column=False):
+    num_columns = len(table_data[0])
+    total_width = 0.9
+    first_column_width_str = ''
+
+    if first_column_width is not None:
+      first_column_width_str = str(first_column_width) + 'cm'
+      total_width -= first_column_width / 16.0
+
+    if caption:
+      self.content.append(r'\begin{table}[htbp]')
+      self.content.append(r'\caption{' + caption + r'}')
+      self.content.append(r'\centering')
+
+    column_format = '|p{' + first_column_width_str + r'}|' + '|'.join(
+      ['X'] * (num_columns - 1)) + '|' if first_column_width is not None else '|' + '|'.join(['X'] * num_columns) + '|'
+    self.content.append(r'\begin{tabularx}{\textwidth}{' + column_format + '}')
+    self.content.append(r'\hline')
+
+    previous_first_column = None
+    rowspan_count = 1
+
+    for i, row in enumerate(table_data):
+      if merge_first_column and i > 0 and row[0] == previous_first_column:
+        rowspan_count += 1
+        self.content.append(' & '.join(['' if j == 0 else cell for j, cell in enumerate(row)]) + r' \\')
+      else:
+        if merge_first_column and i > 0 and rowspan_count > 1:
+          self.content[-rowspan_count] = self.content[-rowspan_count].replace(r'\multirow{1}',
+                                                                              r'\multirow{' + str(rowspan_count) + '}')
+          rowspan_count = 1
+        if merge_first_column and i < len(table_data) - 1 and row[0] == table_data[i + 1][0]:
+          self.content.append(r'\multirow{1}{*}{' + row[0] + '}' + ' & ' + ' & '.join(row[1:]) + r' \\')
+        else:
+          self.content.append(' & '.join(row) + r' \\')
+        previous_first_column = row[0]
+      self.content.append(r'\hline')
+
+    if merge_first_column and rowspan_count > 1:
+      self.content[-rowspan_count] = self.content[-rowspan_count].replace(r'\multirow{1}',
+                                                                          r'\multirow{' + str(rowspan_count) + '}')
+
+    self.content.append(r'\end{tabularx}')
+
+    if caption:
+      self.content.append(r'\end{table}')
+
+  def add_image(self, image_path, caption=None, placement='ht'):
+    if caption:
+      self.content.append(r'\begin{figure}[' + placement + r']')
+      self.content.append(r'\centering')
+      self.content.append(r'\includegraphics[width=\textwidth]{' + image_path + r'}')
+      self.content.append(r'\caption{' + caption + r'}')
+      self.content.append(r'\end{figure}')
+    else:
+      self.content.append(r'\begin{figure}[' + placement + r']')
+      self.content.append(r'\centering')
+      self.content.append(r'\includegraphics[width=\textwidth]{' + image_path + r'}')
+      self.content.append(r'\end{figure}')
+
+  def add_itemize(self, items):
+    self.content.append(r'\begin{itemize}')
+    for item in items:
+      self.content.append(r'\item ' + item)
+    self.content.append(r'\end{itemize}')
+
+  def save_report(self):
+    self.content.append(r'\end{document}')
+    with open(self.file_path, 'w') as f:
+      f.write('\n'.join(self.content))
+
+  def compile_to_pdf(self):
+    subprocess.run(['pdflatex', '-output-directory', str(self.output_path), str(self.file_path)])
+

hub/data/energy_systems/montreal_custom_systems.xml

@@ -22,7 +22,7 @@
         </equipment>
         <equipment id="5" type="cooler" fuel_type="electricity">
             <name>Air cooled DX with external condenser</name>
-            <cooling_efficiency>3.23</cooling_efficiency>
+            <cooling_efficiency>2.5</cooling_efficiency>
             <storage>false</storage>
         </equipment>
         <equipment id="6" type="electricity generator" fuel_type="renewable">
@@ -32,8 +32,8 @@
         </equipment>
         <equipment id="7" type="heat pump" fuel_type="electricity">
             <name>Heat Pump</name>
-            <heating_efficiency>2.79</heating_efficiency>
+            <heating_efficiency>2.5</heating_efficiency>
-            <cooling_efficiency>3.23</cooling_efficiency>
+            <cooling_efficiency>3</cooling_efficiency>
             <storage>false</storage>
         </equipment>
     </generation_equipments>
@@ -198,7 +198,7 @@
             <equipments>
                 <generation_id>3</generation_id>
                 <distribution_id>8</distribution_id>
-g            </equipments>
+           </equipments>
         </system>
         <system id="5">
             <name>Single zone packaged rooftop unit with electrical resistance furnace and baseboards and fuel boiler for acs</name>

hub/data/energy_systems/montreal_future_systems.xml

@@ -911,7 +911,7 @@
             <nominal_cooling_output/>
             <minimum_cooling_output/>
             <maximum_cooling_output/>
-            <cooling_efficiency>4.5</cooling_efficiency>
+            <cooling_efficiency>4</cooling_efficiency>
             <electricity_efficiency/>
             <source_temperature/>
             <source_mass_flow/>
@@ -1411,7 +1411,7 @@
 		</demands>
 		<components>
 			<generation_id>23</generation_id>
-			<generation_id>16</generation_id>
+			<generation_id>17</generation_id>
 		</components>
 	</system>
     <system>

monthly_dhw.py

@@ -0,0 +1,56 @@
+import json
+from pathlib import Path
+import matplotlib.pyplot as plt
+import numpy as np
+from matplotlib.ticker import MaxNLocator
+
+output_path = (Path(__file__).parent / 'out_files').resolve()
+
+# File paths for the three JSON files
+file1 = output_path / 'base_case_buildings_data.json'
+file2 = output_path / 'air_to_air_hp_buildings_data.json'
+file3 = output_path / 'air_to_water_hp_buildings_data.json'
+
+# Opening and reading all three JSON files at the same time
+with open(file1) as f1, open(file2) as f2, open(file3) as f3:
+    base_case = json.load(f1)
+    air = json.load(f2)
+    water = json.load(f3)
+
+month_names = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']
+x = np.arange(len(month_names))  # the label locations
+width = 0.25  # the width of the bars
+
+# Prettier colors for each scenario
+colors = ['#66B2FF', '#e74c3c']  # Blue, Red, Green
+
+
+# Plotting heating data for all buildings in a 2x5 grid
+fig, axes = plt.subplots(2, 5, figsize=(20, 10), dpi=96)
+fig.suptitle('Monthly DHW Consumption Comparison Across Buildings', fontsize=16, weight='bold', alpha=0.8)
+axes = axes.flatten()
+
+for idx, building_name in enumerate(base_case.keys()):
+    heating_data = [list(data["monthly_dhw_consumption_kWh"].values()) for data in
+                    [base_case[building_name], water[building_name]]]
+
+    ax = axes[idx]
+    for i, data in enumerate(heating_data):
+        ax.bar(x + (i - 1) * width, data, width, label=f'Scenario {i+1}', color=colors[i], zorder=2)
+
+    # Grid settings
+    ax.grid(which="major", axis='x', color='#DAD8D7', alpha=0.5, zorder=1)
+    ax.grid(which="major", axis='y', color='#DAD8D7', alpha=0.5, zorder=1)
+
+    # Axis labels and title
+    ax.set_title(building_name, fontsize=14, weight='bold', alpha=0.8, pad=10)
+    ax.set_xticks(x)
+    ax.set_xticklabels(month_names, rotation=45, ha='right')
+    ax.xaxis.set_major_locator(MaxNLocator(integer=True))
+    ax.yaxis.set_major_locator(MaxNLocator(integer=True))
+    if idx % 5 == 0:
+        ax.set_ylabel('DHW Consumption (kWh)', fontsize=12, labelpad=10)
+
+fig.legend(['Base Case', 'Scenario 1&2'], loc='upper right', ncol=3)
+plt.tight_layout(rect=[0, 0.03, 1, 0.95])
+plt.savefig(output_path / 'monthly_dhw.png')

monthly_hvac_plots.py

@@ -0,0 +1,87 @@
+import json
+from pathlib import Path
+import matplotlib.pyplot as plt
+import numpy as np
+from matplotlib.ticker import MaxNLocator
+
+output_path = (Path(__file__).parent / 'out_files').resolve()
+
+# File paths for the three JSON files
+file1 = output_path / 'base_case_buildings_data.json'
+file2 = output_path / 'air_to_air_hp_buildings_data.json'
+file3 = output_path / 'air_to_water_hp_buildings_data.json'
+
+# Opening and reading all three JSON files at the same time
+with open(file1) as f1, open(file2) as f2, open(file3) as f3:
+    base_case = json.load(f1)
+    air = json.load(f2)
+    water = json.load(f3)
+
+month_names = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']
+x = np.arange(len(month_names))  # the label locations
+width = 0.25  # the width of the bars
+
+# Prettier colors for each scenario
+colors = ['#66B2FF', '#e74c3c', '#2ecc71']  # Blue, Red, Green
+
+
+# Plotting heating data for all buildings in a 2x5 grid
+fig, axes = plt.subplots(2, 5, figsize=(20, 10), dpi=96)
+fig.suptitle('Monthly Heating Consumption Comparison Across Buildings', fontsize=16, weight='bold', alpha=0.8)
+axes = axes.flatten()
+
+for idx, building_name in enumerate(base_case.keys()):
+    heating_data = [list(data["monthly_heating_consumption_kWh"].values()) for data in
+                    [base_case[building_name], air[building_name], water[building_name]]]
+
+    ax = axes[idx]
+    for i, data in enumerate(heating_data):
+        ax.bar(x + (i - 1) * width, data, width, label=f'Scenario {i+1}', color=colors[i], zorder=2)
+
+    # Grid settings
+    ax.grid(which="major", axis='x', color='#DAD8D7', alpha=0.5, zorder=1)
+    ax.grid(which="major", axis='y', color='#DAD8D7', alpha=0.5, zorder=1)
+
+    # Axis labels and title
+    ax.set_title(building_name, fontsize=14, weight='bold', alpha=0.8, pad=10)
+    ax.set_xticks(x)
+    ax.set_xticklabels(month_names, rotation=45, ha='right')
+    ax.xaxis.set_major_locator(MaxNLocator(integer=True))
+    ax.yaxis.set_major_locator(MaxNLocator(integer=True))
+    if idx % 5 == 0:
+        ax.set_ylabel('Heating Consumption (kWh)', fontsize=12, labelpad=10)
+
+fig.legend(['Base Case', 'Scenario 1', 'Scenario 2'], loc='upper right', ncol=3)
+plt.tight_layout(rect=[0, 0.03, 1, 0.95])
+plt.savefig(output_path / 'monthly_heating.png')
+
+# Plotting cooling data for all buildings in a 2x5 grid
+# Plotting cooling data for all buildings in a 2x5 grid
+fig, axes = plt.subplots(2, 5, figsize=(20, 10), dpi=96)
+fig.suptitle('Monthly Cooling Consumption Comparison Across Buildings', fontsize=16, weight='bold', alpha=0.8)
+axes = axes.flatten()
+
+for idx, building_name in enumerate(base_case.keys()):
+    cooling_data = [list(data["monthly_cooling_consumption_kWh"].values()) for data in
+                    [base_case[building_name], air[building_name], water[building_name]]]
+
+    ax = axes[idx]
+    for i, data in enumerate(cooling_data):
+        ax.bar(x + (i - 1) * width, data, width, label=f'Scenario {i+1}', color=colors[i], zorder=2)
+
+    # Grid settings
+    ax.grid(which="major", axis='x', color='#DAD8D7', alpha=0.5, zorder=1)
+    ax.grid(which="major", axis='y', color='#DAD8D7', alpha=0.5, zorder=1)
+
+    # Axis labels and title
+    ax.set_title(building_name, fontsize=14, weight='bold', alpha=0.8, pad=10)
+    ax.set_xticks(x)
+    ax.set_xticklabels(month_names, rotation=45, ha='right')
+    ax.xaxis.set_major_locator(MaxNLocator(integer=True))
+    ax.yaxis.set_major_locator(MaxNLocator(integer=True))
+    if idx % 5 == 0:
+        ax.set_ylabel('Cooling Consumption (kWh)', fontsize=12, labelpad=10)
+
+fig.legend(['Base Case', 'Scenario 1', 'Scenario 2'], loc='upper right', ncol=3)
+plt.tight_layout(rect=[0, 0.03, 1, 0.95])
+plt.savefig(output_path / 'monthly_cooling.png')

peak_load.py

@@ -0,0 +1,73 @@
+import json
+from pathlib import Path
+import matplotlib.pyplot as plt
+import numpy as np
+from matplotlib.ticker import MaxNLocator
+from matplotlib.patches import Patch
+
+output_path = (Path(__file__).parent / 'out_files').resolve()
+
+# File paths for the three JSON files
+file1 = output_path / 'base_case_buildings_data.json'
+file2 = output_path / 'air_to_air_hp_buildings_data.json'
+file3 = output_path / 'air_to_water_hp_buildings_data.json'
+
+# Opening and reading all three JSON files at the same time
+with open(file1) as f1, open(file2) as f2, open(file3) as f3:
+    base_case = json.load(f1)
+    air = json.load(f2)
+    water = json.load(f3)
+
+month_names = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']
+x = np.arange(len(month_names))  # the label locations
+
+
+# Scenario labels and color palette
+scenarios = ['Scenario 1', 'Scenario 2']
+colors = ['#66B2FF', '#e74c3c']  # Blue for Scenario 1, Red for Scenario 2
+width = 0.25  # Width for each bar
+
+# Creating the grid for peak load comparisons across buildings
+fig, axes = plt.subplots(2, 5, figsize=(20, 10), dpi=96)
+fig.suptitle('Yearly Heating and Cooling Peak Load Comparison Across Buildings', fontsize=16, weight='bold', alpha=0.8)
+axes = axes.flatten()
+
+for idx, building_name in enumerate(base_case.keys()):
+    # Extracting heating and cooling peak loads for each scenario
+    heating_peak_load = [
+        air[building_name]["heating_peak_load_kW"],
+        water[building_name]["heating_peak_load_kW"]
+    ]
+    cooling_peak_load = [
+        air[building_name]["cooling_peak_load_kW"],
+        water[building_name]["cooling_peak_load_kW"]
+    ]
+
+    ax = axes[idx]
+    x = np.arange(2)  # X locations for the "Heating" and "Cooling" groups
+
+    # Plotting each scenario for heating and cooling
+    for i in range(len(scenarios)):
+        ax.bar(x[0] - width + i * width, heating_peak_load[i], width, color=colors[i], zorder=2)
+        ax.bar(x[1] - width + i * width, cooling_peak_load[i], width, color=colors[i], zorder=2)
+
+    # Grid and styling
+    ax.grid(which="major", axis='x', color='#DAD8D7', alpha=0.5, zorder=1)
+    ax.grid(which="major", axis='y', color='#DAD8D7', alpha=0.5, zorder=1)
+
+    # Axis and title settings
+    ax.set_title(building_name, fontsize=14, weight='bold', alpha=0.8, pad=10)
+    ax.set_xticks(x)
+    ax.set_xticklabels(['Heating Peak Load', 'Cooling Peak Load'])
+    ax.xaxis.set_major_locator(MaxNLocator(integer=True))
+    ax.yaxis.set_major_locator(MaxNLocator(integer=True))
+    if idx % 5 == 0:
+        ax.set_ylabel('Peak Load (kW)', fontsize=12, labelpad=10)
+
+# Custom legend handles to ensure color match with scenarios
+legend_handles = [Patch(color=colors[i], label=scenarios[i]) for i in range(len(scenarios))]
+
+# Global legend and layout adjustments
+fig.legend(handles=legend_handles, loc='upper right', ncol=1)
+plt.tight_layout(rect=[0, 0.03, 1, 0.95])
+plt.savefig(output_path / 'peak_loads.png')

scripts/energy_system_sizing_and_simulation_factory.py

@@ -0,0 +1,35 @@
+"""
+EnergySystemSizingSimulationFactory retrieve the energy system archetype sizing and simulation module
+SPDX - License - Identifier: LGPL - 3.0 - or -later
+Copyright © 2022 Concordia CERC group
+Project Coder Saeed Ranjbar saeed.ranjbar@mail.concordia.ca
+"""
+
+from scripts.system_simulation_models.archetype13 import Archetype13
+
+
+
+class EnergySystemsSimulationFactory:
+  """
+  EnergySystemsFactory class
+  """
+
+  def __init__(self, handler, building, output_path):
+    self._output_path = output_path
+    self._handler = '_' + handler.lower()
+    self._building = building
+
+  def _archetype13(self):
+    """
+    Enrich the city by using the sizing and simulation model developed for archetype13 of montreal_future_systems
+    """
+    Archetype13(self._building, self._output_path).enrich_buildings()
+    self._building.level_of_detail.energy_systems = 2
+    self._building.level_of_detail.energy_systems = 2
+
+  def enrich(self):
+    """
+    Enrich the city given to the class using the class given handler
+    :return: None
+    """
+    getattr(self, self._handler, lambda: None)()

scripts/system_simulation_models/archetype13.py

@@ -0,0 +1,391 @@
+import math
+import hub.helpers.constants as cte
+import csv
+from hub.helpers.monthly_values import MonthlyValues
+
+
+class Archetype13:
+  def __init__(self, building, output_path):
+    self._building = building
+    self._name = building.name
+    self._pv_system = building.energy_systems[0]
+    self._hvac_system = building.energy_systems[1]
+    self._dhw_system = building.energy_systems[-1]
+    self._dhw_peak_flow_rate = (building.thermal_zones_from_internal_zones[0].total_floor_area *
+                                building.thermal_zones_from_internal_zones[0].domestic_hot_water.peak_flow *
+                                cte.WATER_DENSITY)
+    self._heating_peak_load = building.heating_peak_load[cte.YEAR][0]
+    self._cooling_peak_load = building.cooling_peak_load[cte.YEAR][0]
+    self._domestic_hot_water_peak_load = building.domestic_hot_water_peak_load[cte.YEAR][0]
+    self._hourly_heating_demand = [demand / cte.HOUR_TO_SECONDS for demand in building.heating_demand[cte.HOUR]]
+    self._hourly_cooling_demand = [demand / cte.HOUR_TO_SECONDS for demand in building.cooling_demand[cte.HOUR]]
+    self._hourly_dhw_demand = [demand / cte.WATTS_HOUR_TO_JULES for demand in
+                               building.domestic_hot_water_heat_demand[cte.HOUR]]
+    self._output_path = output_path
+    self._t_out = building.external_temperature[cte.HOUR]
+    self.results = {}
+    self.dt = 900
+
+  def hvac_sizing(self):
+    storage_factor = 1.5
+    heat_pump = self._hvac_system.generation_systems[1]
+    boiler = self._hvac_system.generation_systems[0]
+    thermal_storage = boiler.energy_storage_systems[0]
+    heat_pump.nominal_heat_output = round(self._heating_peak_load)
+    heat_pump.nominal_cooling_output = round(self._cooling_peak_load)
+    boiler.nominal_heat_output = 0
+    thermal_storage.volume = round(
+      (self._heating_peak_load * storage_factor * cte.WATTS_HOUR_TO_JULES) /
+      (cte.WATER_HEAT_CAPACITY * cte.WATER_DENSITY * 25))
+    return heat_pump, boiler, thermal_storage
+
+  def dhw_sizing(self):
+    storage_factor = 3
+    dhw_hp = self._dhw_system.generation_systems[0]
+    dhw_hp.nominal_heat_output = 0.7 * self._domestic_hot_water_peak_load
+    dhw_hp.source_temperature = self._t_out
+    dhw_tes = dhw_hp.energy_storage_systems[0]
+    dhw_tes.volume = round(
+      (self._domestic_hot_water_peak_load * storage_factor * 3600) / (cte.WATER_HEAT_CAPACITY * cte.WATER_DENSITY * 10))
+    if dhw_tes.volume == 0:
+      dhw_tes.volume = 1
+    return dhw_hp, dhw_tes
+
+  def heating_system_simulation(self):
+    hp, boiler, tes = self.hvac_sizing()
+    heat_efficiency = float(hp.heat_efficiency)
+    cop_curve_coefficients = [float(coefficient) for coefficient in hp.heat_efficiency_curve.coefficients]
+    number_of_ts = int(cte.HOUR_TO_SECONDS / self.dt)
+    demand = [0] + [x for x in self._hourly_heating_demand for _ in range(number_of_ts)]
+    t_out = [0] + [x for x in self._t_out for _ in range(number_of_ts)]
+    hp.source_temperature = self._t_out
+    variable_names = ["t_sup_hp", "t_tank", "t_ret", "m_ch", "m_dis", "q_hp", "q_boiler", "hp_cop",
+                      "hp_electricity", "boiler_gas_consumption", "t_sup_boiler", "boiler_energy_consumption",
+                      "heating_consumption"]
+    num_hours = len(demand)
+    variables = {name: [0] * num_hours for name in variable_names}
+    (t_sup_hp, t_tank, t_ret, m_ch, m_dis, q_hp, q_boiler, hp_cop,
+     hp_electricity, boiler_gas_consumption, t_sup_boiler, boiler_energy_consumption, heating_consumption) = \
+        [variables[name] for name in variable_names]
+    t_tank[0] = 55
+    hp_heating_cap = hp.nominal_heat_output
+    boiler_heating_cap = boiler.nominal_heat_output
+    hp_delta_t = 7
+    boiler_efficiency = float(boiler.heat_efficiency)
+    v, h = float(tes.volume), float(tes.height)
+    r_tot = sum(float(layer.thickness) / float(layer.material.conductivity) for layer in
+                tes.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)
+    # storage temperature prediction
+    for i in range(len(demand) - 1):
+      t_tank[i + 1] = (t_tank[i] +
+                       (m_ch[i] * (t_sup_boiler[i] - t_tank[i]) +
+                        (ua * (t_out[i] - t_tank[i])) / cte.WATER_HEAT_CAPACITY -
+                        m_dis[i] * (t_tank[i] - t_ret[i])) * (self.dt / (cte.WATER_DENSITY * v)))
+      # hp operation
+      if t_tank[i + 1] < 40:
+        q_hp[i + 1] = hp_heating_cap
+        m_ch[i + 1] = q_hp[i + 1] / (cte.WATER_HEAT_CAPACITY * hp_delta_t)
+        t_sup_hp[i + 1] = (q_hp[i + 1] / (m_ch[i + 1] * cte.WATER_HEAT_CAPACITY)) + t_tank[i + 1]
+      elif 40 <= t_tank[i + 1] < 55 and q_hp[i] == 0:
+        q_hp[i + 1] = 0
+        m_ch[i + 1] = 0
+        t_sup_hp[i + 1] = t_tank[i + 1]
+      elif 40 <= t_tank[i + 1] < 55 and q_hp[i] > 0:
+        q_hp[i + 1] = hp_heating_cap
+        m_ch[i + 1] = q_hp[i + 1] / (cte.WATER_HEAT_CAPACITY * hp_delta_t)
+        t_sup_hp[i + 1] = (q_hp[i + 1] / (m_ch[i + 1] * cte.WATER_HEAT_CAPACITY)) + t_tank[i + 1]
+      else:
+        q_hp[i + 1], m_ch[i + 1], t_sup_hp[i + 1] = 0, 0, t_tank[i + 1]
+      t_sup_hp_fahrenheit = 1.8 * t_sup_hp[i + 1] + 32
+      t_out_fahrenheit = 1.8 * t_out[i + 1] + 32
+      if q_hp[i + 1] > 0:
+        hp_cop[i + 1] = (cop_curve_coefficients[0] +
+                         cop_curve_coefficients[1] * t_sup_hp_fahrenheit +
+                         cop_curve_coefficients[2] * t_sup_hp_fahrenheit ** 2 +
+                         cop_curve_coefficients[3] * t_out_fahrenheit +
+                         cop_curve_coefficients[4] * t_out_fahrenheit ** 2 +
+                         cop_curve_coefficients[5] * t_sup_hp_fahrenheit * t_out_fahrenheit)
+        hp_electricity[i + 1] = q_hp[i + 1] / heat_efficiency
+      else:
+        hp_cop[i + 1] = 0
+        hp_electricity[i + 1] = 0
+    # boiler operation
+      if q_hp[i + 1] > 0:
+        if t_sup_hp[i + 1] < 45:
+          q_boiler[i + 1] = boiler_heating_cap
+        elif demand[i + 1] > 0.5 * self._heating_peak_load / self.dt:
+          q_boiler[i + 1] = 0.5 * boiler_heating_cap
+        boiler_energy_consumption[i + 1] = q_boiler[i + 1] / boiler_efficiency
+        # boiler_gas_consumption[i + 1] = (q_boiler[i + 1] * self.dt) / (boiler_efficiency * cte.NATURAL_GAS_LHV)
+        t_sup_boiler[i + 1] = t_sup_hp[i + 1] + (q_boiler[i + 1] / (m_ch[i + 1] * cte.WATER_HEAT_CAPACITY))
+    # storage discharging
+      if demand[i + 1] == 0:
+        m_dis[i + 1] = 0
+        t_ret[i + 1] = t_tank[i + 1]
+      else:
+        if demand[i + 1] > 0.5 * self._heating_peak_load / cte.HOUR_TO_SECONDS:
+          factor = 8
+        else:
+          factor = 4
+        m_dis[i + 1] = self._heating_peak_load / (cte.WATER_HEAT_CAPACITY * factor * cte.HOUR_TO_SECONDS)
+        t_ret[i + 1] = t_tank[i + 1] - demand[i + 1] / (m_dis[i + 1] * cte.WATER_HEAT_CAPACITY)
+    tes.temperature = []
+    hp_electricity_j = [(x * cte.WATTS_HOUR_TO_JULES) / number_of_ts for x in hp_electricity]
+    boiler_consumption_j = [(x * cte.WATTS_HOUR_TO_JULES) / number_of_ts for x in boiler_energy_consumption]
+    hp_hourly = []
+    boiler_hourly = []
+    boiler_sum = 0
+    hp_sum = 0
+    for i in range(1, len(demand)):
+      hp_sum += hp_electricity_j[i]
+      boiler_sum += boiler_consumption_j[i]
+      if (i - 1) % number_of_ts == 0:
+        tes.temperature.append(t_tank[i])
+        hp_hourly.append(hp_sum)
+        boiler_hourly.append(boiler_sum)
+        hp_sum = 0
+        boiler_sum = 0
+    hp.energy_consumption[cte.HEATING] = {}
+    hp.energy_consumption[cte.HEATING][cte.HOUR] = hp_hourly
+    hp.energy_consumption[cte.HEATING][cte.MONTH] = MonthlyValues.get_total_month(
+      hp.energy_consumption[cte.HEATING][cte.HOUR])
+    hp.energy_consumption[cte.HEATING][cte.YEAR] = [
+      sum(hp.energy_consumption[cte.HEATING][cte.MONTH])]
+    boiler.energy_consumption[cte.HEATING] = {}
+    boiler.energy_consumption[cte.HEATING][cte.HOUR] = boiler_hourly
+    boiler.energy_consumption[cte.HEATING][cte.MONTH] = MonthlyValues.get_total_month(
+      boiler.energy_consumption[cte.HEATING][cte.HOUR])
+    boiler.energy_consumption[cte.HEATING][cte.YEAR] = [
+      sum(boiler.energy_consumption[cte.HEATING][cte.MONTH])]
+
+    self.results['Heating Demand (W)'] = demand
+    self.results['HP Heat Output (W)'] = q_hp
+    self.results['HP Source Temperature'] = t_out
+    self.results['HP Supply Temperature'] = t_sup_hp
+    self.results['HP COP'] = hp_cop
+    self.results['HP Electricity Consumption (W)'] = hp_electricity
+    self.results['Boiler Heat Output (W)'] = q_boiler
+    self.results['Boiler Supply Temperature'] = t_sup_boiler
+    self.results['Boiler Gas Consumption'] = boiler_gas_consumption
+    self.results['TES Temperature'] = t_tank
+    self.results['TES Charging Flow Rate (kg/s)'] = m_ch
+    self.results['TES Discharge Flow Rate (kg/s)'] = m_dis
+    self.results['Heating Loop Return Temperature'] = t_ret
+    return hp_hourly, boiler_hourly
+
+  def cooling_system_simulation(self):
+    hp = self.hvac_sizing()[0]
+    eer_curve_coefficients = [float(coefficient) for coefficient in hp.cooling_efficiency_curve.coefficients]
+    cooling_efficiency = float(hp.cooling_efficiency)
+    number_of_ts = int(cte.HOUR_TO_SECONDS / self.dt)
+    demand = [0] + [x for x in self._hourly_cooling_demand for _ in range(number_of_ts)]
+    t_out = [0] + [x for x in self._t_out for _ in range(number_of_ts)]
+    hp.source_temperature = self._t_out
+    variable_names = ["t_sup_hp", "t_ret", "m", "q_hp", "hp_electricity", "hp_eer"]
+    num_hours = len(demand)
+    variables = {name: [0] * num_hours for name in variable_names}
+    (t_sup_hp, t_ret, m, q_hp, hp_electricity, hp_eer) = [variables[name] for name in variable_names]
+    t_ret[0] = 13
+
+    for i in range(1, len(demand)):
+      if demand[i] > 0:
+        m[i] = self._cooling_peak_load / (cte.WATER_HEAT_CAPACITY * 5 * cte.HOUR_TO_SECONDS)
+        if t_ret[i - 1] >= 13:
+          if demand[i] < 0.25 * self._cooling_peak_load / cte.HOUR_TO_SECONDS:
+            q_hp[i] = 0.25 * hp.nominal_cooling_output
+          elif demand[i] < 0.5 * self._cooling_peak_load / cte.HOUR_TO_SECONDS:
+            q_hp[i] = 0.5 * hp.nominal_cooling_output
+          else:
+            q_hp[i] = hp.nominal_cooling_output
+          t_sup_hp[i] = t_ret[i - 1] - q_hp[i] / (m[i] * cte.WATER_HEAT_CAPACITY)
+        else:
+          q_hp[i] = 0
+          t_sup_hp[i] = t_ret[i - 1]
+        if m[i] == 0:
+          t_ret[i] = t_sup_hp[i]
+        else:
+          t_ret[i] = t_sup_hp[i] + demand[i] / (m[i] * cte.WATER_HEAT_CAPACITY)
+      else:
+        m[i] = 0
+        q_hp[i] = 0
+        t_sup_hp[i] = t_ret[i -1]
+        t_ret[i] = t_ret[i - 1]
+      t_sup_hp_fahrenheit = 1.8 * t_sup_hp[i] + 32
+      t_out_fahrenheit = 1.8 * t_out[i] + 32
+      if q_hp[i] > 0:
+        hp_eer[i] = (eer_curve_coefficients[0] +
+                     eer_curve_coefficients[1] * t_sup_hp_fahrenheit +
+                     eer_curve_coefficients[2] * t_sup_hp_fahrenheit ** 2 +
+                     eer_curve_coefficients[3] * t_out_fahrenheit +
+                     eer_curve_coefficients[4] * t_out_fahrenheit ** 2 +
+                     eer_curve_coefficients[5] * t_sup_hp_fahrenheit * t_out_fahrenheit)
+        hp_electricity[i] = q_hp[i] / cooling_efficiency
+      else:
+        hp_eer[i] = 0
+        hp_electricity[i] = 0
+    hp_electricity_j = [(x * cte.WATTS_HOUR_TO_JULES) / number_of_ts for x in hp_electricity]
+    hp_hourly = []
+    hp_sum = 0
+    for i in range(1, len(demand)):
+      hp_sum += hp_electricity_j[i]
+      if (i - 1) % number_of_ts == 0:
+        hp_hourly.append(hp_sum)
+        hp_sum = 0
+    hp.energy_consumption[cte.COOLING] = {}
+    hp.energy_consumption[cte.COOLING][cte.HOUR] = hp_hourly
+    hp.energy_consumption[cte.COOLING][cte.MONTH] = MonthlyValues.get_total_month(
+      hp.energy_consumption[cte.COOLING][cte.HOUR])
+    hp.energy_consumption[cte.COOLING][cte.YEAR] = [
+      sum(hp.energy_consumption[cte.COOLING][cte.MONTH])]
+    self.results['Cooling Demand (W)'] = demand
+    self.results['HP Cooling Output (W)'] = q_hp
+    self.results['HP Cooling Supply Temperature'] = t_sup_hp
+    self.results['HP Cooling COP'] = hp_eer
+    self.results['HP Electricity Consumption'] = hp_electricity
+    self.results['Cooling Loop Flow Rate (kg/s)'] = m
+    self.results['Cooling Loop Return Temperature'] = t_ret
+    return hp_hourly
+
+  def dhw_system_simulation(self):
+    hp, tes = self.dhw_sizing()
+    heat_efficiency = float(hp.heat_efficiency)
+    cop_curve_coefficients = [float(coefficient) for coefficient in hp.heat_efficiency_curve.coefficients]
+    number_of_ts = int(cte.HOUR_TO_SECONDS / self.dt)
+    demand = [0] + [x for x in self._hourly_dhw_demand for _ in range(number_of_ts)]
+    t_out = [0] + [x for x in self._t_out for _ in range(number_of_ts)]
+    variable_names = ["t_sup_hp", "t_tank", "m_ch", "m_dis", "q_hp", "q_coil", "hp_cop",
+                      "hp_electricity", "available hot water (m3)", "refill flow rate (kg/s)"]
+    num_hours = len(demand)
+    variables = {name: [0] * num_hours for name in variable_names}
+    (t_sup_hp, t_tank, m_ch, m_dis, m_refill, q_hp, q_coil, hp_cop, hp_electricity, v_dhw) = \
+        [variables[name] for name in variable_names]
+    t_tank[0] = 70
+    v_dhw[0] = tes.volume
+
+    hp_heating_cap = hp.nominal_heat_output
+    hp_delta_t = 8
+    v, h = float(tes.volume), float(tes.height)
+    r_tot = sum(float(layer.thickness) / float(layer.material.conductivity) for layer in
+                tes.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)
+    freshwater_temperature = 18
+    for i in range(len(demand) - 1):
+      delta_t_demand = demand[i] * (self.dt / (cte.WATER_DENSITY * cte.WATER_HEAT_CAPACITY * v))
+      if t_tank[i] < 65:
+        q_hp[i] = hp_heating_cap
+      delta_t_hp = q_hp[i] * (self.dt / (cte.WATER_DENSITY * cte.WATER_HEAT_CAPACITY * v))
+      if demand[i] > 0:
+        dhw_needed = (demand[i] * cte.HOUR_TO_SECONDS) / (cte.WATER_HEAT_CAPACITY * t_tank[i] * cte.WATER_DENSITY)
+        m_dis[i] = dhw_needed * cte.WATER_DENSITY / cte.HOUR_TO_SECONDS
+        m_refill[i] = m_dis[i]
+      delta_t_freshwater = m_refill[i] * (t_tank[i] - freshwater_temperature) * (self.dt / (v * cte.WATER_DENSITY))
+      diff = delta_t_freshwater + delta_t_demand - delta_t_hp
+      if diff > 0:
+          if diff > 0:
+            power = diff * (cte.WATER_DENSITY * cte.WATER_HEAT_CAPACITY * v) / self.dt
+            if power <= float(tes.heating_coil_capacity):
+              q_coil[i] = power
+            else:
+              q_coil[i] = float(tes.heating_coil_capacity)
+      delta_t_coil = q_coil[i] * (self.dt / (cte.WATER_DENSITY * cte.WATER_HEAT_CAPACITY * v))
+
+      if q_hp[i] > 0:
+        m_ch[i] = q_hp[i] / (cte.WATER_HEAT_CAPACITY * hp_delta_t)
+        t_sup_hp[i] = (q_hp[i] / (m_ch[i] * cte.WATER_HEAT_CAPACITY)) + t_tank[i]
+      else:
+        m_ch[i] = 0
+        t_sup_hp[i] = t_tank[i]
+      t_sup_hp_fahrenheit = 1.8 * t_sup_hp[i] + 32
+      t_out_fahrenheit = 1.8 * t_out[i] + 32
+      if q_hp[i] > 0:
+        hp_cop[i] = (cop_curve_coefficients[0] +
+                     cop_curve_coefficients[1] * t_sup_hp_fahrenheit +
+                     cop_curve_coefficients[2] * t_sup_hp_fahrenheit ** 2 +
+                     cop_curve_coefficients[3] * t_out_fahrenheit +
+                     cop_curve_coefficients[4] * t_out_fahrenheit ** 2 +
+                     cop_curve_coefficients[5] * t_sup_hp_fahrenheit * t_out_fahrenheit)
+        hp_electricity[i] = q_hp[i] / heat_efficiency
+      else:
+        hp_cop[i] = 0
+        hp_electricity[i] = 0
+
+      t_tank[i + 1] = t_tank[i] + (delta_t_hp - delta_t_freshwater - delta_t_demand + delta_t_coil)
+    tes.temperature = []
+    hp_electricity_j = [(x * cte.WATTS_HOUR_TO_JULES) / number_of_ts for x in hp_electricity]
+    heating_coil_j = [(x * cte.WATTS_HOUR_TO_JULES) / number_of_ts for x in q_coil]
+    hp_hourly = []
+    coil_hourly = []
+    coil_sum = 0
+    hp_sum = 0
+    for i in range(1, len(demand)):
+      hp_sum += hp_electricity_j[i]
+      coil_sum += heating_coil_j[i]
+      if (i - 1) % number_of_ts == 0:
+        tes.temperature.append(t_tank[i])
+        hp_hourly.append(hp_sum)
+        coil_hourly.append(coil_sum)
+        hp_sum = 0
+        coil_sum = 0
+
+    hp.energy_consumption[cte.DOMESTIC_HOT_WATER] = {}
+    hp.energy_consumption[cte.DOMESTIC_HOT_WATER][cte.HOUR] = hp_hourly
+    hp.energy_consumption[cte.DOMESTIC_HOT_WATER][cte.MONTH] = MonthlyValues.get_total_month(
+      hp.energy_consumption[cte.DOMESTIC_HOT_WATER][cte.HOUR])
+    hp.energy_consumption[cte.DOMESTIC_HOT_WATER][cte.YEAR] = [
+      sum(hp.energy_consumption[cte.DOMESTIC_HOT_WATER][cte.MONTH])]
+    tes.heating_coil_energy_consumption = {}
+    tes.heating_coil_energy_consumption[cte.HOUR] = coil_hourly
+    tes.heating_coil_energy_consumption[cte.MONTH] = MonthlyValues.get_total_month(
+      tes.heating_coil_energy_consumption[cte.HOUR])
+    tes.heating_coil_energy_consumption[cte.YEAR] = [
+      sum(tes.heating_coil_energy_consumption[cte.MONTH])]
+    tes.temperature = t_tank
+
+    self.results['DHW Demand (W)'] = demand
+    self.results['DHW HP Heat Output (W)'] = q_hp
+    self.results['DHW HP Electricity Consumption (W)'] = hp_electricity
+    self.results['DHW HP Source Temperature'] = t_out
+    self.results['DHW HP Supply Temperature'] = t_sup_hp
+    self.results['DHW HP COP'] = hp_cop
+    self.results['DHW TES Heating Coil Heat Output (W)'] = q_coil
+    self.results['DHW TES Temperature'] = t_tank
+    self.results['DHW TES Charging Flow Rate (kg/s)'] = m_ch
+    self.results['DHW Flow Rate (kg/s)'] = m_dis
+    self.results['DHW TES Refill Flow Rate (kg/s)'] = m_refill
+    self.results['Available Water in Tank (m3)'] = v_dhw
+    return hp_hourly, coil_hourly
+
+  def enrich_buildings(self):
+    hp_heating, boiler_consumption = self.heating_system_simulation()
+    hp_cooling = self.cooling_system_simulation()
+    hp_dhw, heating_coil = self.dhw_system_simulation()
+    heating_consumption = [hp_heating[i] + boiler_consumption[i] for i in range(len(hp_heating))]
+    dhw_consumption = [hp_dhw[i] + heating_coil[i] for i in range(len(hp_dhw))]
+    self._building.heating_consumption[cte.HOUR] = 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])]
+    self._building.cooling_consumption[cte.HOUR] = hp_cooling
+    self._building.cooling_consumption[cte.MONTH] = (
+      MonthlyValues.get_total_month(self._building.cooling_consumption[cte.HOUR]))
+    self._building.cooling_consumption[cte.YEAR] = [sum(self._building.cooling_consumption[cte.MONTH])]
+    self._building.domestic_hot_water_consumption[cte.HOUR] = dhw_consumption
+    self._building.domestic_hot_water_consumption[cte.MONTH] = (
+      MonthlyValues.get_total_month(self._building.domestic_hot_water_consumption[cte.HOUR]))
+    self._building.domestic_hot_water_consumption[cte.YEAR] = [sum(self._building.domestic_hot_water_consumption[cte.MONTH])]
+    file_name = f'energy_system_simulation_results_{self._name}.csv'
+    with open(self._output_path / file_name, 'w', newline='') as csvfile:
+      output_file = csv.writer(csvfile)
+      # Write header
+      output_file.writerow(self.results.keys())
+      # Write data
+      output_file.writerows(zip(*self.results.values()))