feature: add district heating network creator

Reviewed-on: https://nextgenerations-cities.encs.concordia.ca/gitea/s_ranjbar/energy_system_modelling_workflow/pulls/16

.idea/energy_system_modelling_workflow.iml

@@ -2,7 +2,7 @@
 <module type="PYTHON_MODULE" version="4">
   <component name="NewModuleRootManager">
     <content url="file://$MODULE_DIR$" />
-    <orderEntry type="jdk" jdkName="hub" jdkType="Python SDK" />
+    <orderEntry type="jdk" jdkName="C:\Users\marie.proville\miniconda3\envs\Hub" jdkType="Python SDK" />
     <orderEntry type="sourceFolder" forTests="false" />
   </component>
 </module>

.idea/misc.xml

@@ -1,4 +1,7 @@
 <?xml version="1.0" encoding="UTF-8"?>
 <project version="4">
+  <component name="Black">
+    <option name="sdkUUID" value="97386509-ec9b-4dd7-929b-7585219c0447" />
+  </component>
   <component name="ProjectRootManager" version="2" project-jdk-name="hub" project-jdk-type="Python SDK" />
 </project>

energy_system_retrofit.py

@@ -1,4 +1,3 @@
-from scripts.geojson_creator import process_geojson
 from pathlib import Path
 import subprocess
 from scripts.ep_run_enrich import energy_plus_workflow
@@ -8,57 +7,94 @@ from hub.imports.construction_factory import ConstructionFactory
 from hub.imports.usage_factory import UsageFactory
 from hub.imports.weather_factory import WeatherFactory
 from hub.imports.results_factory import ResultFactory
-from scripts.energy_system_analysis_report import EnergySystemAnalysisReport
+from scripts.energy_system_retrofit_report import EnergySystemRetrofitReport
+from scripts.geojson_creator import process_geojson
 from scripts import random_assignation
 from hub.imports.energy_systems_factory import EnergySystemsFactory
 from scripts.energy_system_sizing import SystemSizing
-from scripts.energy_system_retrofit_results import system_results, new_system_results
+from scripts.solar_angles import CitySolarAngles
+from scripts.pv_sizing_and_simulation import PVSizingSimulation
+from scripts.energy_system_retrofit_results import consumption_data, cost_data
 from scripts.energy_system_sizing_and_simulation_factory import EnergySystemsSimulationFactory
 from scripts.costs.cost import Cost
-from scripts.costs.constants import SKIN_RETROFIT_AND_SYSTEM_RETROFIT_AND_PV, SYSTEM_RETROFIT_AND_PV
+from scripts.costs.constants import SKIN_RETROFIT_AND_SYSTEM_RETROFIT_AND_PV, SYSTEM_RETROFIT_AND_PV, CURRENT_STATUS
 import hub.helpers.constants as cte
 from hub.exports.exports_factory import ExportsFactory
+from scripts.pv_feasibility import pv_feasibility
+
 # Specify the GeoJSON file path
+input_files_path = (Path(__file__).parent / 'input_files')
+input_files_path.mkdir(parents=True, exist_ok=True)
 geojson_file = process_geojson(x=-73.5681295982132, y=45.49218262677643, diff=0.0001)
-file_path = (Path(__file__).parent / 'input_files' / 'output_buildings.geojson')
-# Specify the output path for the PDF file
+geojson_file_path = input_files_path / 'output_buildings.geojson'
 output_path = (Path(__file__).parent / 'out_files').resolve()
-# Create city object from GeoJSON file
-city = GeometryFactory('geojson',
-                       path=file_path,
+output_path.mkdir(parents=True, exist_ok=True)
+energy_plus_output_path = output_path / 'energy_plus_outputs'
+energy_plus_output_path.mkdir(parents=True, exist_ok=True)
+simulation_results_path = (Path(__file__).parent / 'out_files' / 'simulation_results').resolve()
+simulation_results_path.mkdir(parents=True, exist_ok=True)
+sra_output_path = output_path / 'sra_outputs'
+sra_output_path.mkdir(parents=True, exist_ok=True)
+cost_analysis_output_path = output_path / 'cost_analysis'
+cost_analysis_output_path.mkdir(parents=True, exist_ok=True)
+city = GeometryFactory(file_type='geojson',
+                       path=geojson_file_path,
                        height_field='height',
                        year_of_construction_field='year_of_construction',
                        function_field='function',
                        function_to_hub=Dictionaries().montreal_function_to_hub_function).city
-# Enrich city data
 ConstructionFactory('nrcan', city).enrich()
-
 UsageFactory('nrcan', city).enrich()
 WeatherFactory('epw', city).enrich()
-ExportsFactory('sra', city, output_path).export()
-sra_path = (output_path / f'{city.name}_sra.xml').resolve()
+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, output_path).enrich()
-energy_plus_workflow(city)
+ResultFactory('sra', city, sra_output_path).enrich()
+pv_feasibility(-73.5681295982132, 45.49218262677643, 0.0001, selected_buildings=city.buildings)
+energy_plus_workflow(city, energy_plus_output_path)
+solar_angles = CitySolarAngles(city.name,
+                               city.latitude,
+                               city.longitude,
+                               tilt_angle=45,
+                               surface_azimuth_angle=180).calculate
 random_assignation.call_random(city.buildings, random_assignation.residential_systems_percentage)
 EnergySystemsFactory('montreal_custom', city).enrich()
 SystemSizing(city.buildings).montreal_custom()
-current_system = new_system_results(city.buildings)
+current_status_energy_consumption = consumption_data(city)
+current_status_life_cycle_cost = {}
+for building in city.buildings:
+  cost_retrofit_scenario = CURRENT_STATUS
+  lcc_dataframe = Cost(building=building,
+                       retrofit_scenario=cost_retrofit_scenario,
+                       fuel_tariffs=['Electricity-D', 'Gas-Energir']).life_cycle
+  lcc_dataframe.to_csv(cost_analysis_output_path / f'{building.name}_current_status_lcc.csv')
+  current_status_life_cycle_cost[f'{building.name}'] = cost_data(building, lcc_dataframe, cost_retrofit_scenario)
 random_assignation.call_random(city.buildings, random_assignation.residential_new_systems_percentage)
 EnergySystemsFactory('montreal_future', city).enrich()
 for building in city.buildings:
-  EnergySystemsSimulationFactory('archetype1', building=building, output_path=output_path).enrich()
-  print(building.energy_consumption_breakdown[cte.ELECTRICITY][cte.COOLING] +
-        building.energy_consumption_breakdown[cte.ELECTRICITY][cte.HEATING] +
-        building.energy_consumption_breakdown[cte.ELECTRICITY][cte.DOMESTIC_HOT_WATER])
-new_system = new_system_results(city.buildings)
-# EnergySystemAnalysisReport(city, output_path).create_report(current_system, new_system)
+  if 'PV' in building.energy_systems_archetype_name:
+    ghi = [x / cte.WATTS_HOUR_TO_JULES for x in building.roofs[0].global_irradiance[cte.HOUR]]
+    pv_sizing_simulation = PVSizingSimulation(building,
+                                              solar_angles,
+                                              tilt_angle=45,
+                                              module_height=1,
+                                              module_width=2,
+                                              ghi=ghi)
+    pv_sizing_simulation.pv_output()
+  if building.energy_systems_archetype_name == 'PV+4Pipe+DHW':
+    EnergySystemsSimulationFactory('archetype13', building=building, output_path=simulation_results_path).enrich()
+retrofitted_energy_consumption = consumption_data(city)
+retrofitted_life_cycle_cost = {}
 for building in city.buildings:
-  costs = Cost(building=building, retrofit_scenario=SYSTEM_RETROFIT_AND_PV).life_cycle
-  costs.to_csv(output_path / f'{building.name}_lcc.csv')
-  (costs.loc['global_operational_costs', f'Scenario {SYSTEM_RETROFIT_AND_PV}'].
-   to_csv(output_path / f'{building.name}_op.csv'))
-  costs.loc['global_capital_costs', f'Scenario {SYSTEM_RETROFIT_AND_PV}'].to_csv(
-    output_path / f'{building.name}_cc.csv')
-  costs.loc['global_maintenance_costs', f'Scenario {SYSTEM_RETROFIT_AND_PV}'].to_csv(
-    output_path / f'{building.name}_m.csv')
+  cost_retrofit_scenario = SYSTEM_RETROFIT_AND_PV
+  lcc_dataframe = Cost(building=building,
+                       retrofit_scenario=cost_retrofit_scenario,
+                       fuel_tariffs=['Electricity-D', 'Gas-Energir']).life_cycle
+  lcc_dataframe.to_csv(cost_analysis_output_path / f'{building.name}_retrofitted_lcc.csv')
+  retrofitted_life_cycle_cost[f'{building.name}'] = cost_data(building, lcc_dataframe, cost_retrofit_scenario)
+for i in range(12):
+  dhw_consumption = 0
+  for building in city.buildings:
+    dhw_consumption += building.domestic_hot_water_consumption[cte.MONTH][i] / 3.6e6
+
+

hub/catalog_factories/cost/montreal_complete_cost_catalog.py

@@ -48,7 +48,7 @@ class MontrealNewCatalog(Catalog):
       construction = float(archetype['incomes']['subsidies']['construction']['#text'])
       hvac = float(archetype['incomes']['subsidies']['hvac']['#text'])
       photovoltaic_system = float(archetype['incomes']['subsidies']['photovoltaic']['#text'])
-      electricity_exports = float(archetype['incomes']['electricity_export']['#text']) / 1000 / 3600
+      electricity_exports = float(archetype['incomes']['electricity_export']['#text'])
       reduction_tax = float(archetype['incomes']['tax_reduction']['#text']) / 100
       income = Income(construction_subsidy=construction,
                       hvac_subsidy=hvac,

hub/catalog_factories/data_models/energy_systems/emission_system.py

@@ -10,7 +10,7 @@ class EmissionSystem:
   """
   Emission system class
   """
-  def __init__(self, system_id, model_name=None, system_type=None, parasitic_energy_consumption=None):
+  def __init__(self, system_id, model_name=None, system_type=None, parasitic_energy_consumption=0):
 
     self._system_id = system_id
     self._model_name = model_name

hub/catalog_factories/energy_systems/montreal_custom_catalog.py

@@ -135,7 +135,7 @@ class MontrealCustomCatalog(Catalog):
       equipment_id = float(equipment['@id'])
       equipment_type = equipment['@type']
       model_name = equipment['name']
-      parasitic_consumption = None
+      parasitic_consumption = 0
       if 'parasitic_consumption' in equipment:
         parasitic_consumption = float(equipment['parasitic_consumption']['#text']) / 100
 

hub/catalog_factories/energy_systems/montreal_future_system_catalogue.py

@@ -262,7 +262,7 @@ class MontrealFutureSystemCatalogue(Catalog):
         system_id = None
         model_name = None
         system_type = None
-        parasitic_energy_consumption = None
+        parasitic_energy_consumption = 0
         emission_system = EmissionSystem(system_id=system_id,
                                          model_name=model_name,
                                          system_type=system_type,

hub/city_model_structure/building.py

@@ -484,7 +484,7 @@ class Building(CityObject):
       monthly_values = PeakLoads().peak_loads_from_hourly(self.domestic_hot_water_heat_demand[cte.HOUR])
     if monthly_values is None:
       return None
-    results[cte.MONTH] = [x for x in monthly_values]
+    results[cte.MONTH] = [x / cte.WATTS_HOUR_TO_JULES for x in monthly_values]
     results[cte.YEAR] = [max(monthly_values) / cte.WATTS_HOUR_TO_JULES]
     return results
 
@@ -810,39 +810,16 @@ class Building(CityObject):
     Get total electricity produced onsite in J
     return: dict
     """
-    orientation_losses_factor = {cte.MONTH: {'north': [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
-                                             'east': [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
-                                             'south': [2.137931, 1.645503, 1.320946, 1.107817, 0.993213, 0.945175,
-                                                       0.967949, 1.065534, 1.24183, 1.486486, 1.918033, 2.210526],
-                                             'west': [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]},
-                                 cte.YEAR: {'north': [0],
-                                            'east': [0],
-                                            'south': [1.212544],
-                                            'west': [0]}
-                                 }
-
-    # Add other systems whenever new ones appear
-    if self.energy_systems is None:
-      return self._onsite_electrical_production
-    for energy_system in self.energy_systems:
-      for generation_system in energy_system.generation_systems:
-        if generation_system.system_type == cte.PHOTOVOLTAIC:
-          if generation_system.electricity_efficiency is not None:
-            _efficiency = float(generation_system.electricity_efficiency)
-          else:
-            _efficiency = 0
-          self._onsite_electrical_production = {}
-          for _key in self.roofs[0].global_irradiance.keys():
-            _results = [0 for _ in range(0, len(self.roofs[0].global_irradiance[_key]))]
-            for surface in self.roofs:
-              if _key in orientation_losses_factor:
-                _results = [x + y * _efficiency * surface.perimeter_area
-                            * surface.solar_collectors_area_reduction_factor * z
-                            for x, y, z in zip(_results, surface.global_irradiance[_key],
-                                               orientation_losses_factor[_key]['south'])]
-            self._onsite_electrical_production[_key] = _results
     return self._onsite_electrical_production
 
+  @onsite_electrical_production.setter
+  def onsite_electrical_production(self, value):
+    """
+    set onsite electrical production from external pv simulations
+    :return:
+    """
+    self._onsite_electrical_production = value
+
   @property
   def lower_corner(self):
     """
@@ -913,10 +890,3 @@ class Building(CityObject):
     self._fuel_consumption_breakdown = fuel_breakdown
     return self._fuel_consumption_breakdown
 
-  @property
-  def pv_generation(self) -> dict:
-    return self._pv_generation
-
-  @pv_generation.setter
-  def pv_generation(self, value):
-    self._pv_generation = value

hub/city_model_structure/energy_systems/emission_system.py

@@ -13,7 +13,7 @@ class EmissionSystem:
   def __init__(self):
     self._model_name = None
     self._type = None
-    self._parasitic_energy_consumption = None
+    self._parasitic_energy_consumption = 0
 
   @property
   def model_name(self):

hub/data/costs/montreal_costs_completed.xml

@@ -187,7 +187,7 @@
 				<hvac cost_unit="%">1.5</hvac>
 				<photovoltaic cost_unit="%">3.6</photovoltaic>
 			</subsidies>
-            <electricity_export cost_unit="currency/kWh">0.07</electricity_export>
+            <electricity_export cost_unit="currency/kWh">0.075</electricity_export>
             <tax_reduction cost_unit="%">5</tax_reduction>
         </incomes>
     </archetype>

hub/data/energy_systems/montreal_future_systems.xml

@@ -927,7 +927,7 @@
 			    <dependant_variable>COP</dependant_variable>
 		        <parameters>source_temperature</parameters>
 		        <parameters>supply_temperature</parameters>
-		        <coefficients  a="-0.000277" b="0.019639" c="0.000004" d="0.012190" e="-0.00010" f="-0.000277"/>
+		        <coefficients  a="1.039924" b="0.014600" c="0.000006" d="-0.05026" e="0.000635" f="-0.000154"/>
             </heat_efficiency_curve>
             <cooling_output_curve/>
             <cooling_fuel_consumption_curve/>
@@ -1052,7 +1052,7 @@
             <nominal_heat_output/>
             <minimum_heat_output/>
             <maximum_heat_output/>
-            <heat_efficiency>3.5</heat_efficiency>
+            <heat_efficiency>3.2</heat_efficiency>
             <reversible/>
             <fuel_type>electricity</fuel_type>
             <source_medium>Air</source_medium>
@@ -1076,7 +1076,7 @@
 			    <dependant_variable>COP</dependant_variable>
 		        <parameters>source_temperature</parameters>
 		        <parameters>supply_temperature</parameters>
-		        <coefficients  a="1.990668" b="0" c="0" d="-0.027252" e="0.000131" f="0"/>
+		        <coefficients  a="1.19713" b="0.077849" c="-0.0000016" d="-0.02675" e="0.000296" f="-0.00112"/>
             </heat_efficiency_curve>
             <cooling_output_curve/>
             <cooling_fuel_consumption_curve/>
@@ -1438,6 +1438,29 @@
 			<generation_id>27</generation_id>
 		</components>
 	</system>
+    <system>
+		<id>11</id>
+		<name>Central Heating System َASHP Gas-Boiler TES</name>
+		<schema>schemas/ASHP+TES+GasBoiler.jpg</schema>
+		<demands>
+			<demand>heating</demand>
+		</demands>
+		<components>
+			<generation_id>23</generation_id>
+			<generation_id>16</generation_id>
+		</components>
+	</system>
+    <system>
+		<id>12</id>
+		<name>Unitary ASHP Cooling System</name>
+		<schema>schemas/ASHP+TES+GasBoiler.jpg</schema>
+		<demands>
+			<demand>cooling</demand>
+		</demands>
+		<components>
+			<generation_id>23</generation_id>
+		</components>
+	</system>
   </systems>
 
   <system_archetypes>
@@ -1528,6 +1551,23 @@
             <system_id>10</system_id>
 		</systems>
 	</system_archetype>
+	<system_archetype id="14">
+		<name>Central Heating+Unitary Cooling+Unitary DHW</name>
+		<systems>
+			<system_id>10</system_id>
+			<system_id>11</system_id>
+            <system_id>12</system_id>
+		</systems>
+	</system_archetype>
+	<system_archetype id="15">
+		<name>Central Heating+Unitary Cooling+Unitary DHW+PV</name>
+		<systems>
+			<system_id>7</system_id>
+            <system_id>10</system_id>
+			<system_id>11</system_id>
+            <system_id>12</system_id>
+		</systems>
+	</system_archetype>
   </system_archetypes>
 </EnergySystemCatalog>
 

hub/imports/construction/nrcan_physics_parameters.py

@@ -32,10 +32,21 @@ class NrcanPhysicsParameters:
     city = self._city
     nrcan_catalog = ConstructionCatalogFactory('nrcan').catalog
     for building in city.buildings:
-      if building.function not in Dictionaries().hub_function_to_nrcan_construction_function:
-        logging.error('Building %s has an unknown building function %s', building.name, building.function)
+      main_function = None
+      functions = building.function.split('_')
+      if len(functions) > 1:
+        maximum_percentage = 0
+        for function in functions:
+          percentage_and_function = function.split('-')
+          if float(percentage_and_function[0]) > maximum_percentage:
+            maximum_percentage = float(percentage_and_function[0])
+            main_function = percentage_and_function[-1]
+      else:
+        main_function = functions[-1]
+      if main_function not in Dictionaries().hub_function_to_nrcan_construction_function:
+        logging.error('Building %s has an unknown building function %s', building.name, main_function)
         continue
-      function = Dictionaries().hub_function_to_nrcan_construction_function[building.function]
+      function = Dictionaries().hub_function_to_nrcan_construction_function[main_function]
       try:
         archetype = self._search_archetype(nrcan_catalog, function, building.year_of_construction, self._climate_zone)
 

hub/imports/energy_systems/montreal_custom_energy_system_parameters.py

@@ -136,10 +136,14 @@ class MontrealCustomEnergySystemParameters:
       _distribution_system.distribution_consumption_variable_flow = \
         archetype_distribution_system.distribution_consumption_variable_flow
       _distribution_system.heat_losses = archetype_distribution_system.heat_losses
-      _emission_system = None
+      _generic_emission_system = None
       if archetype_distribution_system.emission_systems is not None:
-        _emission_system = EmissionSystem()
-      _distribution_system.emission_systems = [_emission_system]
+        _emission_systems = []
+        for emission_system in archetype_distribution_system.emission_systems:
+          _generic_emission_system = EmissionSystem()
+          _generic_emission_system.parasitic_energy_consumption = emission_system.parasitic_energy_consumption
+          _emission_systems.append(_generic_emission_system)
+        _distribution_system.emission_systems = _emission_systems
       _distribution_systems.append(_distribution_system)
     return _distribution_systems
 

hub/imports/energy_systems/montreal_future_energy_systems_parameters.py

@@ -185,10 +185,14 @@ class MontrealFutureEnergySystemParameters:
         _distribution_system.distribution_consumption_variable_flow = \
           archetype_distribution_system.distribution_consumption_variable_flow
         _distribution_system.heat_losses = archetype_distribution_system.heat_losses
-        _emission_system = None
+        _generic_emission_system = None
         if archetype_distribution_system.emission_systems is not None:
-          _emission_system = EmissionSystem()
-        _distribution_system.emission_systems = [_emission_system]
+          _emission_systems = []
+          for emission_system in archetype_distribution_system.emission_systems:
+            _generic_emission_system = EmissionSystem()
+            _generic_emission_system.parasitic_energy_consumption = emission_system.parasitic_energy_consumption
+            _emission_systems.append(_generic_emission_system)
+          _distribution_system.emission_systems = _emission_systems
         _distribution_systems.append(_distribution_system)
       return _distribution_systems
 

hub/imports/geometry/geojson.py

@@ -127,6 +127,27 @@ class Geojson:
         function = None
         if self._function_field is not None:
           function = str(feature['properties'][self._function_field])
+          if function == 'Mixed use' or function == 'mixed use':
+            function_parts = []
+            if 'usages' in feature['properties']:
+              usages = feature['properties']['usages']
+              for usage in usages:
+                if self._function_to_hub is not None and usage['usage'] in self._function_to_hub:
+                  function_parts.append(f"{usage['percentage']}-{self._function_to_hub[usage['usage']]}")
+                else:
+                  function_parts.append(f"{usage['percentage']}-{usage['usage']}")
+            else:
+              for key, value in feature['properties'].items():
+                if key.startswith("mixed_type_") and not key.endswith("_percentage"):
+                  type_key = key
+                  percentage_key = f"{key}_percentage"
+                  if percentage_key in feature['properties']:
+                    if self._function_to_hub is not None and feature['properties'][type_key] in self._function_to_hub:
+                      usage_function = self._function_to_hub[feature['properties'][type_key]]
+                      function_parts.append(f"{feature['properties'][percentage_key]}-{usage_function}")
+                    else:
+                      function_parts.append(f"{feature['properties'][percentage_key]}-{feature['properties'][type_key]}")
+            function = "_".join(function_parts)
           if self._function_to_hub is not None:
             # use the transformation dictionary to retrieve the proper function
             if function in self._function_to_hub:

hub/imports/usage/comnet_usage_parameters.py

@@ -35,29 +35,60 @@ class ComnetUsageParameters:
     city = self._city
     comnet_catalog = UsageCatalogFactory('comnet').catalog
     for building in city.buildings:
-      usage_name = Dictionaries().hub_usage_to_comnet_usage[building.function]
-      try:
-        archetype_usage = self._search_archetypes(comnet_catalog, usage_name)
-      except KeyError:
-        logging.error('Building %s has unknown usage archetype for usage %s', building.name, usage_name)
-        continue
-
-      for internal_zone in building.internal_zones:
-        if internal_zone.area is None:
-          raise TypeError('Internal zone area not defined, ACH cannot be calculated')
-        if internal_zone.volume is None:
-          raise TypeError('Internal zone volume not defined, ACH cannot be calculated')
-        if internal_zone.area <= 0:
-          raise TypeError('Internal zone area is zero, ACH cannot be calculated')
-        volume_per_area = internal_zone.volume / internal_zone.area
-        usage = Usage()
-        usage.name = usage_name
-        self._assign_values(usage, archetype_usage, volume_per_area, building.cold_water_temperature)
-        usage.percentage = 1
-        self._calculate_reduced_values_from_extended_library(usage, archetype_usage)
-
-        internal_zone.usages = [usage]
+      usages = []
+      comnet_archetype_usages = []
+      building_functions = building.function.split('_')
+      for function in building_functions:
+        usages.append(function.split('-'))
+      for usage in usages:
+        comnet_usage_name = Dictionaries().hub_usage_to_comnet_usage[usage[-1]]
+        try:
+          comnet_archetype_usage = self._search_archetypes(comnet_catalog, comnet_usage_name)
+          comnet_archetype_usages.append(comnet_archetype_usage)
+        except KeyError:
+          logging.error('Building %s has unknown usage archetype for usage %s', building.name, comnet_usage_name)
+          continue
+      for (i, internal_zone) in enumerate(building.internal_zones):
+        internal_zone_usages = []
+        if len(building.internal_zones) > 1:
+          volume_per_area = 0
+          if internal_zone.area is None:
+            logging.error('Building %s has internal zone area not defined, ACH cannot be calculated for usage %s',
+                          building.name, usages[i][-1])
+            continue
+          if internal_zone.volume is None:
+            logging.error('Building %s has internal zone volume not defined, ACH cannot be calculated for usage %s',
+                          building.name, usages[i][-1])
+            continue
+          if internal_zone.area <= 0:
+            logging.error('Building %s has internal zone area equal to 0, ACH cannot be calculated for usage %s',
+                          building.name,  usages[i][-1])
+            continue
+          volume_per_area += internal_zone.volume / internal_zone.area
+          usage = Usage()
+          usage.name = usages[i][-1]
+          self._assign_values(usage, comnet_archetype_usages[i], volume_per_area, building.cold_water_temperature)
+          usage.percentage = 1
+          self._calculate_reduced_values_from_extended_library(usage, comnet_archetype_usages[i])
+          internal_zone_usages.append(usage)
+        else:
+          if building.storeys_above_ground is None:
+            logging.error('Building %s no number of storeys assigned, ACH cannot be calculated for usage %s',
+                          building.name,  usages)
+            continue
+          volume_per_area = building.volume / building.floor_area / building.storeys_above_ground
+          for (j, mixed_usage) in enumerate(usages):
+            usage = Usage()
+            usage.name = mixed_usage[-1]
+            if len(usages) > 1:
+              usage.percentage = float(mixed_usage[0]) / 100
+            else:
+              usage.percentage = 1
+            self._assign_values(usage, comnet_archetype_usages[j], volume_per_area, building.cold_water_temperature)
+            self._calculate_reduced_values_from_extended_library(usage, comnet_archetype_usages[j])
+            internal_zone_usages.append(usage)
 
+        internal_zone.usages = internal_zone_usages
   @staticmethod
   def _search_archetypes(comnet_catalog, usage_name):
     comnet_archetypes = comnet_catalog.entries('archetypes').usages

hub/imports/usage/nrcan_usage_parameters.py

@@ -33,53 +33,72 @@ class NrcanUsageParameters:
     city = self._city
     nrcan_catalog = UsageCatalogFactory('nrcan').catalog
     comnet_catalog = UsageCatalogFactory('comnet').catalog
-
     for building in city.buildings:
-      usage_name = Dictionaries().hub_usage_to_nrcan_usage[building.function]
-      try:
-        archetype_usage = self._search_archetypes(nrcan_catalog, usage_name)
-      except KeyError:
-        logging.error('Building %s has unknown usage archetype for usage %s', building.name, usage_name)
-        continue
+      usages = []
+      nrcan_archetype_usages = []
+      comnet_archetype_usages = []
+      building_functions = building.function.split('_')
+      for function in building_functions:
+        usages.append(function.split('-'))
+      for usage in usages:
+        usage_name = Dictionaries().hub_usage_to_nrcan_usage[usage[-1]]
+        try:
+          archetype_usage = self._search_archetypes(nrcan_catalog, usage_name)
+          nrcan_archetype_usages.append(archetype_usage)
+        except KeyError:
+          logging.error('Building %s has unknown usage archetype for usage %s', building.name, usage_name)
+          continue
+        comnet_usage_name = Dictionaries().hub_usage_to_comnet_usage[usage[-1]]
+        try:
+          comnet_archetype_usage = self._search_archetypes(comnet_catalog, comnet_usage_name)
+          comnet_archetype_usages.append(comnet_archetype_usage)
+        except KeyError:
+          logging.error('Building %s has unknown usage archetype for usage %s', building.name, comnet_usage_name)
+          continue
 
-      comnet_usage_name = Dictionaries().hub_usage_to_comnet_usage[building.function]
-      try:
-        comnet_archetype_usage = self._search_archetypes(comnet_catalog, comnet_usage_name)
-      except KeyError:
-        logging.error('Building %s has unknown usage archetype for usage %s', building.name, comnet_usage_name)
-        continue
-
-      for internal_zone in building.internal_zones:
+      for (i, internal_zone) in enumerate(building.internal_zones):
+        internal_zone_usages = []
         if len(building.internal_zones) > 1:
           volume_per_area = 0
           if internal_zone.area is None:
             logging.error('Building %s has internal zone area not defined, ACH cannot be calculated for usage %s',
-                          building.name, usage_name)
+                          building.name, usages[i][-1])
             continue
           if internal_zone.volume is None:
             logging.error('Building %s has internal zone volume not defined, ACH cannot be calculated for usage %s',
-                          building.name, usage_name)
+                          building.name, usages[i][-1])
             continue
           if internal_zone.area <= 0:
             logging.error('Building %s has internal zone area equal to 0, ACH cannot be calculated for usage %s',
-                          building.name,  usage_name)
+                          building.name,  usages[i][-1])
             continue
           volume_per_area += internal_zone.volume / internal_zone.area
+          usage = Usage()
+          usage.name = usages[i][-1]
+          self._assign_values(usage, nrcan_archetype_usages[i], volume_per_area, building.cold_water_temperature)
+          self._assign_comnet_extra_values(usage, comnet_archetype_usages[i], nrcan_archetype_usages[i].occupancy.occupancy_density)
+          usage.percentage = 1
+          self._calculate_reduced_values_from_extended_library(usage, nrcan_archetype_usages[i])
+          internal_zone_usages.append(usage)
         else:
           if building.storeys_above_ground is None:
             logging.error('Building %s no number of storeys assigned, ACH cannot be calculated for usage %s',
-                          building.name,  usage_name)
+                          building.name,  usages)
             continue
           volume_per_area = building.volume / building.floor_area / building.storeys_above_ground
+          for (j, mixed_usage) in enumerate(usages):
+            usage = Usage()
+            usage.name = mixed_usage[-1]
+            if len(usages) > 1:
+              usage.percentage = float(mixed_usage[0]) / 100
+            else:
+              usage.percentage = 1
+            self._assign_values(usage, nrcan_archetype_usages[j], volume_per_area, building.cold_water_temperature)
+            self._assign_comnet_extra_values(usage, comnet_archetype_usages[j], nrcan_archetype_usages[j].occupancy.occupancy_density)
+            self._calculate_reduced_values_from_extended_library(usage, nrcan_archetype_usages[j])
+            internal_zone_usages.append(usage)
 
-        usage = Usage()
-        usage.name = usage_name
-        self._assign_values(usage, archetype_usage, volume_per_area, building.cold_water_temperature)
-        self._assign_comnet_extra_values(usage, comnet_archetype_usage, archetype_usage.occupancy.occupancy_density)
-        usage.percentage = 1
-        self._calculate_reduced_values_from_extended_library(usage, archetype_usage)
-
-        internal_zone.usages = [usage]
+        internal_zone.usages = internal_zone_usages
 
   @staticmethod
   def _search_archetypes(catalog, usage_name):

input_files/output_buildings_expanded.geojson

@@ -0,0 +1,863 @@
+{
+  "type": "FeatureCollection",
+  "features": [
+    {
+      "type": "Feature",
+      "geometry": {
+        "type": "Polygon",
+        "coordinates": [
+          [
+            [
+              -73.56769087843276,
+              45.49251875903776
+            ],
+            [
+              -73.56765050367694,
+              45.492560280202284
+            ],
+            [
+              -73.5677794213865,
+              45.49262188364245
+            ],
+            [
+              -73.56781916241786,
+              45.49258006136105
+            ],
+            [
+              -73.56769087843276,
+              45.49251875903776
+            ]
+          ]
+        ]
+      },
+      "id": 173347,
+      "properties": {
+        "name": "01044617",
+        "address": "rue Victor-Hugo (MTL) 1666",
+        "function": "1000",
+        "height": 9,
+        "year_of_construction": 1986
+      }
+    },
+    {
+      "type": "Feature",
+      "geometry": {
+        "type": "Polygon",
+        "coordinates": [
+          [
+            [
+              -73.56765050367694,
+              45.492560280202284
+            ],
+            [
+              -73.56761436875776,
+              45.49259744179384
+            ],
+            [
+              -73.5676075694645,
+              45.49260454199484
+            ],
+            [
+              -73.56773226889548,
+              45.49266394156485
+            ],
+            [
+              -73.56773726906921,
+              45.49266624130272
+            ],
+            [
+              -73.5677794213865,
+              45.49262188364245
+            ],
+            [
+              -73.56765050367694,
+              45.492560280202284
+            ]
+          ]
+        ]
+      },
+      "id": 173348,
+      "properties": {
+        "name": "01044619",
+        "address": "rue Victor-Hugo (MTL) 1670",
+        "function": "1000",
+        "height": 9,
+        "year_of_construction": 1986
+      }
+    },
+    {
+      "type": "Feature",
+      "geometry": {
+        "type": "Polygon",
+        "coordinates": [
+          [
+            [
+              -73.56829026835214,
+              45.492524742569145
+            ],
+            [
+              -73.56849646900322,
+              45.49262354174874
+            ],
+            [
+              -73.56861067001111,
+              45.492505541343576
+            ],
+            [
+              -73.56864076915663,
+              45.492519941474434
+            ],
+            [
+              -73.56866246900178,
+              45.49249754209202
+            ],
+            [
+              -73.56867696946317,
+              45.49250454136644
+            ],
+            [
+              -73.56867726964143,
+              45.49250414255471
+            ],
+            [
+              -73.56881486931461,
+              45.492362042624144
+            ],
+            [
+              -73.56881686903772,
+              45.492359941181455
+            ],
+            [
+              -73.5688004699483,
+              45.49235084193039
+            ],
+            [
+              -73.56882097012145,
+              45.4923320417195
+            ],
+            [
+              -73.56879846891101,
+              45.49232034109352
+            ],
+            [
+              -73.56883736970825,
+              45.492284841271946
+            ],
+            [
+              -73.56886806888434,
+              45.492256240993704
+            ],
+            [
+              -73.56885337003277,
+              45.49224914198001
+            ],
+            [
+              -73.56890226932418,
+              45.49219894164121
+            ],
+            [
+              -73.56851866897392,
+              45.49201434154299
+            ],
+            [
+              -73.56837326884313,
+              45.492163841620254
+            ],
+            [
+              -73.56864696910176,
+              45.49229554163243
+            ],
+            [
+              -73.5685268682051,
+              45.49241904187041
+            ],
+            [
+              -73.56825396962694,
+              45.49228824183907
+            ],
+            [
+              -73.56810906858335,
+              45.49243794104013
+            ],
+            [
+              -73.56829026835214,
+              45.492524742569145
+            ]
+          ]
+        ]
+      },
+      "id": 173403,
+      "properties": {
+        "name": "01044334",
+        "address": "rue Saint-Jacques (MTL) 1460",
+        "function": "1000",
+        "height": 15,
+        "year_of_construction": 1985
+      }
+    },
+    {
+      "type": "Feature",
+      "geometry": {
+        "type": "Polygon",
+        "coordinates": [
+          [
+            [
+              -73.5683896684674,
+              45.491800342137736
+            ],
+            [
+              -73.56838616878639,
+              45.49180414157881
+            ],
+            [
+              -73.56850686988925,
+              45.49185994152571
+            ],
+            [
+              -73.56851286844197,
+              45.4918626410622
+            ],
+            [
+              -73.56855549071014,
+              45.49181750806087
+            ],
+            [
+              -73.56842962331187,
+              45.49175738300567
+            ],
+            [
+              -73.5683896684674,
+              45.491800342137736
+            ]
+          ]
+        ]
+      },
+      "id": 174898,
+      "properties": {
+        "name": "01044590",
+        "address": "rue Victor-Hugo (MTL) 1600",
+        "function": "1000",
+        "height": 9,
+        "year_of_construction": 1986
+      }
+    },
+    {
+      "type": "Feature",
+      "geometry": {
+        "type": "Polygon",
+        "coordinates": [
+          [
+            [
+              -73.5680637695714,
+              45.49212884162544
+            ],
+            [
+              -73.56802228176146,
+              45.49217205619571
+            ],
+            [
+              -73.56815668696326,
+              45.49223626189717
+            ],
+            [
+              -73.56815766959974,
+              45.49223524178655
+            ],
+            [
+              -73.56818746886172,
+              45.49224944155107
+            ],
+            [
+              -73.56822816806918,
+              45.49220694186927
+            ],
+            [
+              -73.5680637695714,
+              45.49212884162544
+            ]
+          ]
+        ]
+      },
+      "id": 175785,
+      "properties": {
+        "name": "01044602",
+        "address": "rue Victor-Hugo (MTL) 1630",
+        "function": "1000",
+        "height": 12,
+        "year_of_construction": 1986
+      }
+    },
+    {
+      "type": "Feature",
+      "geometry": {
+        "type": "Polygon",
+        "coordinates": [
+          [
+            [
+              -73.56850793693103,
+              45.49167318076048
+            ],
+            [
+              -73.56846877951091,
+              45.4917152818903
+            ],
+            [
+              -73.56859506290321,
+              45.491775605518725
+            ],
+            [
+              -73.56863463503653,
+              45.491733702062774
+            ],
+            [
+              -73.56850793693103,
+              45.49167318076048
+            ]
+          ]
+        ]
+      },
+      "id": 175910,
+      "properties": {
+        "name": "01044586",
+        "address": "rue Victor-Hugo (MTL) 1590",
+        "function": "1000",
+        "height": 9,
+        "year_of_construction": 1986
+      }
+    },
+    {
+      "type": "Feature",
+      "geometry": {
+        "type": "Polygon",
+        "coordinates": [
+          [
+            [
+              -73.56817543449134,
+              45.49201384773851
+            ],
+            [
+              -73.56813497596143,
+              45.49205532773507
+            ],
+            [
+              -73.56826745951075,
+              45.492118613912375
+            ],
+            [
+              -73.56830763251781,
+              45.49207699906335
+            ],
+            [
+              -73.56817543449134,
+              45.49201384773851
+            ]
+          ]
+        ]
+      },
+      "id": 176056,
+      "properties": {
+        "name": "01044599",
+        "address": "rue Victor-Hugo (MTL) 1620",
+        "function": "1000",
+        "height": 8,
+        "year_of_construction": 1986
+      }
+    },
+    {
+      "type": "Feature",
+      "geometry": {
+        "type": "Polygon",
+        "coordinates": [
+          [
+            [
+              -73.56772876855176,
+              45.49247194194522
+            ],
+            [
+              -73.56773406949068,
+              45.492474341387755
+            ],
+            [
+              -73.56773125185198,
+              45.492477239659124
+            ],
+            [
+              -73.56785890467093,
+              45.492538239964624
+            ],
+            [
+              -73.56789966910456,
+              45.49249534173201
+            ],
+            [
+              -73.56776616865103,
+              45.49243264153464
+            ],
+            [
+              -73.56772876855176,
+              45.49247194194522
+            ]
+          ]
+        ]
+      },
+      "id": 176261,
+      "properties": {
+        "name": "01044613",
+        "address": "rue Victor-Hugo (MTL) 1656",
+        "function": "1000",
+        "height": 10,
+        "year_of_construction": 1986
+      }
+    },
+    {
+      "type": "Feature",
+      "geometry": {
+        "type": "Polygon",
+        "coordinates": [
+          [
+            [
+              -73.56802228176146,
+              45.49217205619571
+            ],
+            [
+              -73.56798225825526,
+              45.492213743742184
+            ],
+            [
+              -73.56811660206223,
+              45.49227791893211
+            ],
+            [
+              -73.56815668696326,
+              45.49223626189717
+            ],
+            [
+              -73.56802228176146,
+              45.49217205619571
+            ]
+          ]
+        ]
+      },
+      "id": 176293,
+      "properties": {
+        "name": "01044604",
+        "address": "rue Victor-Hugo (MTL) 1636",
+        "function": "1000",
+        "height": 12,
+        "year_of_construction": 1986
+      }
+    },
+    {
+      "type": "Feature",
+      "geometry": {
+        "type": "Polygon",
+        "coordinates": [
+          [
+            [
+              -73.56790222258577,
+              45.49229712328457
+            ],
+            [
+              -73.56785996900595,
+              45.49234104192853
+            ],
+            [
+              -73.56799446861396,
+              45.49240484193282
+            ],
+            [
+              -73.56803643080562,
+              45.49236123475947
+            ],
+            [
+              -73.56790222258577,
+              45.49229712328457
+            ]
+          ]
+        ]
+      },
+      "id": 176296,
+      "properties": {
+        "name": "01044611",
+        "address": "rue Victor-Hugo (MTL) 1650",
+        "function": "1000",
+        "height": 10,
+        "year_of_construction": 1986
+      }
+    },
+    {
+      "type": "Feature",
+      "geometry": {
+        "type": "Polygon",
+        "coordinates": [
+          [
+            [
+              -73.56798225825526,
+              45.492213743742184
+            ],
+            [
+              -73.56794223597048,
+              45.4922554321734
+            ],
+            [
+              -73.56807651582375,
+              45.49231957685336
+            ],
+            [
+              -73.56811660206223,
+              45.49227791893211
+            ],
+            [
+              -73.56798225825526,
+              45.492213743742184
+            ]
+          ]
+        ]
+      },
+      "id": 176298,
+      "properties": {
+        "name": "01044607",
+        "address": "rue Victor-Hugo (MTL) 1640",
+        "function": "1000",
+        "height": 12,
+        "year_of_construction": 1986
+      }
+    },
+    {
+      "type": "Feature",
+      "geometry": {
+        "type": "Polygon",
+        "coordinates": [
+          [
+            [
+              -73.56742736898599,
+              45.49184704208998
+            ],
+            [
+              -73.56761256873325,
+              45.491896142437554
+            ],
+            [
+              -73.56766926915839,
+              45.4917902412014
+            ],
+            [
+              -73.56766956853903,
+              45.49179024192391
+            ],
+            [
+              -73.56792966911675,
+              45.49183254222432
+            ],
+            [
+              -73.56793006788594,
+              45.491831141828406
+            ],
+            [
+              -73.56794526884076,
+              45.49174634219527
+            ],
+            [
+              -73.56794516904765,
+              45.49174634225465
+            ],
+            [
+              -73.56753896905731,
+              45.491638642248425
+            ],
+            [
+              -73.56742736898599,
+              45.49184704208998
+            ]
+          ]
+        ]
+      },
+      "id": 176918,
+      "properties": {
+        "name": "01097185",
+        "address": "rue Victor-Hugo (MTL) 1591",
+        "function": "1000",
+        "height": 10,
+        "year_of_construction": 1987
+      }
+    },
+    {
+      "type": "Feature",
+      "geometry": {
+        "type": "Polygon",
+        "coordinates": [
+          [
+            [
+              -73.56773125185198,
+              45.492477239659124
+            ],
+            [
+              -73.56769087843276,
+              45.49251875903776
+            ],
+            [
+              -73.56781916241786,
+              45.49258006136105
+            ],
+            [
+              -73.56785890467093,
+              45.492538239964624
+            ],
+            [
+              -73.56773125185198,
+              45.492477239659124
+            ]
+          ]
+        ]
+      },
+      "id": 178164,
+      "properties": {
+        "name": "01044615",
+        "address": "rue Victor-Hugo (MTL) 1660",
+        "function": "1000",
+        "height": 9,
+        "year_of_construction": 1986
+      }
+    },
+    {
+      "type": "Feature",
+      "geometry": {
+        "type": "Polygon",
+        "coordinates": [
+          [
+            [
+              -73.56846877951091,
+              45.4917152818903
+            ],
+            [
+              -73.56842962331187,
+              45.49175738300567
+            ],
+            [
+              -73.56855549071014,
+              45.49181750806087
+            ],
+            [
+              -73.56859506290321,
+              45.491775605518725
+            ],
+            [
+              -73.56846877951091,
+              45.4917152818903
+            ]
+          ]
+        ]
+      },
+      "id": 179679,
+      "properties": {
+        "name": "01044588",
+        "address": "rue Victor-Hugo (MTL) 1596",
+        "function": "1000",
+        "height": 9,
+        "year_of_construction": 1986
+      }
+    },
+    {
+      "type": "Feature",
+      "geometry": {
+        "type": "Polygon",
+        "coordinates": [
+          [
+            [
+              -73.56825635009473,
+              45.49193088860213
+            ],
+            [
+              -73.56821589168355,
+              45.491972368627906
+            ],
+            [
+              -73.5683477837006,
+              45.4920353716151
+            ],
+            [
+              -73.56838787594006,
+              45.49199371809223
+            ],
+            [
+              -73.56825635009473,
+              45.49193088860213
+            ]
+          ]
+        ]
+      },
+      "id": 179789,
+      "properties": {
+        "name": "01044595",
+        "address": "rue Victor-Hugo (MTL) 1610",
+        "function": "1000",
+        "height": 8,
+        "year_of_construction": 1986
+      }
+    },
+    {
+      "type": "Feature",
+      "geometry": {
+        "type": "Polygon",
+        "coordinates": [
+          [
+            [
+              -73.56821589168355,
+              45.491972368627906
+            ],
+            [
+              -73.56817543449134,
+              45.49201384773851
+            ],
+            [
+              -73.56830763251781,
+              45.49207699906335
+            ],
+            [
+              -73.5683477837006,
+              45.4920353716151
+            ],
+            [
+              -73.56821589168355,
+              45.491972368627906
+            ]
+          ]
+        ]
+      },
+      "id": 181310,
+      "properties": {
+        "name": "01044597",
+        "address": "rue Victor-Hugo (MTL) 1616",
+        "function": "1000",
+        "height": 8,
+        "year_of_construction": 1986
+      }
+    },
+    {
+      "type": "Feature",
+      "geometry": {
+        "type": "Polygon",
+        "coordinates": [
+          [
+            [
+              -73.56809506939487,
+              45.49209624228538
+            ],
+            [
+              -73.56809246893268,
+              45.4920988416879
+            ],
+            [
+              -73.56821287000538,
+              45.49216124158406
+            ],
+            [
+              -73.56822186852654,
+              45.49216584161625
+            ],
+            [
+              -73.56826745951075,
+              45.492118613912375
+            ],
+            [
+              -73.56813497596143,
+              45.49205532773507
+            ],
+            [
+              -73.56809506939487,
+              45.49209624228538
+            ]
+          ]
+        ]
+      },
+      "id": 182393,
+      "properties": {
+        "name": "01044601",
+        "address": "rue Victor-Hugo (MTL) 1626",
+        "function": "1000",
+        "height": 8,
+        "year_of_construction": 1986
+      }
+    },
+    {
+      "type": "Feature",
+      "geometry": {
+        "type": "Polygon",
+        "coordinates": [
+          [
+            [
+              -73.56790756893894,
+              45.492291541967774
+            ],
+            [
+              -73.56790222258577,
+              45.49229712328457
+            ],
+            [
+              -73.56803643080562,
+              45.49236123475947
+            ],
+            [
+              -73.56807651582375,
+              45.49231957685336
+            ],
+            [
+              -73.56794223597048,
+              45.4922554321734
+            ],
+            [
+              -73.56790756893894,
+              45.492291541967774
+            ]
+          ]
+        ]
+      },
+      "id": 182442,
+      "properties": {
+        "name": "01044609",
+        "address": "rue Victor-Hugo (MTL) 1646",
+        "function": "1000",
+        "height": 11,
+        "year_of_construction": 1986
+      }
+    },
+    {
+      "type": "Feature",
+      "geometry": {
+        "type": "Polygon",
+        "coordinates": [
+          [
+            [
+              -73.56829706912258,
+              45.49188914205178
+            ],
+            [
+              -73.56825635009473,
+              45.49193088860213
+            ],
+            [
+              -73.56838787594006,
+              45.49199371809223
+            ],
+            [
+              -73.56842846901456,
+              45.49195154234486
+            ],
+            [
+              -73.56829706912258,
+              45.49188914205178
+            ]
+          ]
+        ]
+      },
+      "id": 182546,
+      "properties": {
+        "name": "01044592",
+        "address": "rue Victor-Hugo (MTL) 1606",
+        "function": "1000",
+        "height": 8,
+        "year_of_construction": 1986
+      }
+    }
+  ]
+}

input_files/test_geojson.geojson

@@ -0,0 +1,55 @@
+{
+  "type": "FeatureCollection",
+  "features": [
+    {
+      "type": "Feature",
+      "geometry": {
+        "type": "Polygon",
+        "coordinates": [
+          [
+            [
+              -73.58000127109773,
+              45.49613461675315
+            ],
+            [
+              -73.57962787855432,
+              45.496524875557746
+            ],
+            [
+              -73.57996357265695,
+              45.49668114195629
+            ],
+            [
+              -73.57996427397713,
+              45.496680342403664
+            ],
+            [
+              -73.58034707390021,
+              45.49625804233725
+            ],
+            [
+              -73.58034697395713,
+              45.496257942524835
+            ],
+            [
+              -73.58000127109773,
+              45.49613461675315
+            ]
+          ]
+        ]
+      },
+      "id": 179764,
+      "properties": {
+        "name": "01119274",
+        "address": "rue Guy (MTL) 2157",
+        "function": "Mixed use",
+        "mixed_type_1": "commercial",
+        "mixed_type_1_percentage": 50,
+        "mixed_type_2": "6000",
+        "mixed_type_2_percentage": 50,
+        "height": 62,
+        "year_of_construction": 1954
+      }
+    }
+  ]
+}

input_files/test_geojson1.geojson

@@ -0,0 +1,52 @@
+{
+  "type": "FeatureCollection",
+  "features": [
+    {
+      "type": "Feature",
+      "geometry": {
+        "type": "Polygon",
+        "coordinates": [
+          [
+            [
+              -73.58000127109773,
+              45.49613461675315
+            ],
+            [
+              -73.57962787855432,
+              45.496524875557746
+            ],
+            [
+              -73.57996357265695,
+              45.49668114195629
+            ],
+            [
+              -73.57996427397713,
+              45.496680342403664
+            ],
+            [
+              -73.58034707390021,
+              45.49625804233725
+            ],
+            [
+              -73.58034697395713,
+              45.496257942524835
+            ],
+            [
+              -73.58000127109773,
+              45.49613461675315
+            ]
+          ]
+        ]
+      },
+      "id": 179764,
+      "properties": {
+        "name": "01119274",
+        "address": "rue Guy (MTL) 2157",
+        "function": "Mixed use",
+        "usages": [{"usage": "commercial", "percentage": 50}, {"usage": "6000", "percentage": 50}],
+        "height": 62,
+        "year_of_construction": 1954
+      }
+    }
+  ]
+}

main.py

@@ -0,0 +1,86 @@
+from pathlib import Path
+import subprocess
+from scripts.ep_run_enrich import energy_plus_workflow
+from hub.imports.geometry_factory import GeometryFactory
+from hub.helpers.dictionaries import Dictionaries
+from hub.imports.construction_factory import ConstructionFactory
+from hub.imports.usage_factory import UsageFactory
+from hub.imports.weather_factory import WeatherFactory
+from hub.imports.results_factory import ResultFactory
+from scripts.energy_system_retrofit_report import EnergySystemRetrofitReport
+from scripts.geojson_creator import process_geojson
+from scripts import random_assignation
+from hub.imports.energy_systems_factory import EnergySystemsFactory
+from scripts.energy_system_sizing import SystemSizing
+from scripts.solar_angles import CitySolarAngles
+from scripts.pv_sizing_and_simulation import PVSizingSimulation
+from scripts.energy_system_retrofit_results import consumption_data, cost_data
+from scripts.energy_system_sizing_and_simulation_factory import EnergySystemsSimulationFactory
+from scripts.costs.cost import Cost
+from scripts.costs.constants import SKIN_RETROFIT_AND_SYSTEM_RETROFIT_AND_PV, SYSTEM_RETROFIT_AND_PV, CURRENT_STATUS
+import hub.helpers.constants as cte
+from hub.exports.exports_factory import ExportsFactory
+from scripts.pv_feasibility import pv_feasibility
+import matplotlib.pyplot as plt
+import numpy as np
+# Specify the GeoJSON file path
+data = {}
+input_files_path = (Path(__file__).parent / 'input_files')
+input_files_path.mkdir(parents=True, exist_ok=True)
+# geojson_file = process_geojson(x=-73.58001358793511, y=45.496445294438715, diff=0.0001)
+geojson_file_path = input_files_path / 'test_geojson1.geojson'
+output_path = (Path(__file__).parent / 'out_files').resolve()
+output_path.mkdir(parents=True, exist_ok=True)
+energy_plus_output_path = output_path / 'energy_plus_outputs'
+energy_plus_output_path.mkdir(parents=True, exist_ok=True)
+simulation_results_path = (Path(__file__).parent / 'out_files' / 'simulation_results').resolve()
+simulation_results_path.mkdir(parents=True, exist_ok=True)
+sra_output_path = output_path / 'sra_outputs'
+sra_output_path.mkdir(parents=True, exist_ok=True)
+cost_analysis_output_path = output_path / 'cost_analysis'
+cost_analysis_output_path.mkdir(parents=True, exist_ok=True)
+city = GeometryFactory(file_type='geojson',
+                       path=geojson_file_path,
+                       height_field='height',
+                       year_of_construction_field='year_of_construction',
+                       function_field='function',
+                       function_to_hub=Dictionaries().montreal_function_to_hub_function).city
+# ConstructionFactory('nrcan', city).enrich()
+# UsageFactory('nrcan', city).enrich()
+# WeatherFactory('epw', city).enrich()
+# energy_plus_workflow(city, energy_plus_output_path)
+# data[f'{city.buildings[0].function}'] = city.buildings[0].heating_demand[cte.YEAR][0] / 3.6e9
+# city.buildings[0].function = cte.COMMERCIAL
+# ConstructionFactory('nrcan', city).enrich()
+# UsageFactory('nrcan', city).enrich()
+# energy_plus_workflow(city, energy_plus_output_path)
+# data[f'{city.buildings[0].function}'] = city.buildings[0].heating_demand[cte.YEAR][0] / 3.6e9
+# city.buildings[0].function = cte.MEDIUM_OFFICE
+# ConstructionFactory('nrcan', city).enrich()
+# UsageFactory('nrcan', city).enrich()
+# energy_plus_workflow(city, energy_plus_output_path)
+# data[f'{city.buildings[0].function}'] = city.buildings[0].heating_demand[cte.YEAR][0] / 3.6e9
+# categories = list(data.keys())
+# values = list(data.values())
+# # Plotting
+# fig, ax = plt.subplots(figsize=(10, 6), dpi=96)
+# fig.suptitle('Impact of different usages on yearly heating demand', fontsize=16, weight='bold', alpha=.8)
+# ax.bar(categories, values, color=['#2196f3', '#ff5a5f', '#4caf50'], 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('Building Type', fontsize=12, labelpad=10)
+# ax.set_ylabel('Energy Consumption (MWh)', fontsize=14, labelpad=10)
+# ax.yaxis.set_major_locator(plt.MaxNLocator(integer=True))
+# ax.set_xticks(np.arange(len(categories)))
+# ax.set_xticklabels(categories, rotation=45, ha='right')
+# ax.bar_label(ax.containers[0], padding=3, color='black', fontsize=12, rotation=0)
+# ax.spines[['top', 'left', 'bottom']].set_visible(False)
+# ax.spines['right'].set_linewidth(1.1)
+# # Set a white background
+# fig.patch.set_facecolor('white')
+# # Adjust the margins around the plot area
+# plt.subplots_adjust(left=0.1, right=0.9, top=0.85, bottom=0.25)
+# # Save the plot
+# plt.savefig('plot_nrcan.png', bbox_inches='tight')
+# plt.close()
+print('test')

scripts/costs/capital_costs.py

@@ -13,7 +13,7 @@ from hub.city_model_structure.building import Building
 import hub.helpers.constants as cte
 from scripts.costs.configuration import Configuration
 from scripts.costs.constants import (SKIN_RETROFIT, SKIN_RETROFIT_AND_SYSTEM_RETROFIT_AND_PV,
-                                     SYSTEM_RETROFIT_AND_PV, CURRENT_STATUS, PV)
+                                     SYSTEM_RETROFIT_AND_PV, CURRENT_STATUS, PV, SYSTEM_RETROFIT)
 from scripts.costs.cost_base import CostBase
 
 
@@ -32,10 +32,12 @@ class CapitalCosts(CostBase):
         'B3010_opaque_roof',
         'B1010_superstructure',
         'D2010_photovoltaic_system',
-        'D3020_heat_and_cooling_generating_systems',
-        'D3040_distribution_systems',
-        'D3050_other_hvac_ahu',
-        'D3060_storage_systems',
+        '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'
@@ -45,10 +47,12 @@ class CapitalCosts(CostBase):
     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_heat_and_cooling_generating_systems'] = 0
-    self._yearly_capital_costs.loc[0, 'D3040_distribution_systems'] = 0
-    self._yearly_capital_costs.loc[0, 'D3080_other_hvac_ahu'] = 0
-    self._yearly_capital_costs.loc[0, 'D3060_storage_systems'] = 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(
@@ -107,7 +111,7 @@ class CapitalCosts(CostBase):
     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):
+    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
@@ -151,14 +155,20 @@ class CapitalCosts(CostBase):
     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] == 'generation':
+      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]
@@ -170,19 +180,25 @@ class CapitalCosts(CostBase):
 
     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):
-      self._yearly_capital_costs.loc[0, 'D3020_heat_and_cooling_generating_systems'] = (
+    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, 'D3040_distribution_systems'] = (
+      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, 'D3060_storage_systems'] = (
+      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_distribution_equipment,
-            capital_cost_energy_storage_equipment, capital_cost_domestic_hot_water_equipment, capital_cost_lighting, capital_cost_hvac)
+    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')
@@ -202,49 +218,79 @@ class CapitalCosts(CostBase):
           system_investment_costs[0] * self._configuration.percentage_credit
         )
       )
-      self._yearly_capital_costs.loc[year, 'D3020_heat_and_cooling_generating_systems'] = (
+      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, 'D3040_distribution_systems'] = (
+      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, 'D3060_storage_systems'] = (
+      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, 'D40_dhw'] = (
+      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] == 'generation':
-              reposition_cost_heating_and_cooling_equipment = chapter.item(component).reposition[0] * component_sizes[i] * costs_increase
-              self._yearly_capital_costs.loc[year, 'D3020_heat_and_cooling_generating_systems'] += reposition_cost_heating_and_cooling_equipment
+            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
+              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, 'D3040_distribution_systems'] += reposition_cost_distribution_equipment
+              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, 'D3060_storage_systems'] += reposition_cost_energy_storage_equipment
+              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'] += (
@@ -280,8 +326,11 @@ class CapitalCosts(CostBase):
               system_components.append(self.boiler_type(generation_system))
             else:
               system_components.append('D302010_template_heat')
-          elif cte.HEATING or cte.COOLING in demand_types:
-            component_categories.append('generation')
+          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)
@@ -290,11 +339,18 @@ class CapitalCosts(CostBase):
               item_type = self.boiler_type(generation_system)
               system_components.append(item_type)
             else:
-              if cte.COOLING in demand_types and cte.HEATING not in demand_types:
+              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:

scripts/costs/constants.py

@@ -12,10 +12,12 @@ 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
+  PV,
+  SYSTEM_RETROFIT
 ]

scripts/costs/cost.py

@@ -93,10 +93,12 @@ class Cost:
         global_capital_costs['B1010_superstructure']
     )
     df_capital_costs_systems = (
-        global_capital_costs['D3020_heat_and_cooling_generating_systems'] +
-        global_capital_costs['D3040_distribution_systems'] +
-        global_capital_costs['D3050_other_hvac_ahu'] +
-        global_capital_costs['D3060_storage_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']
     )

scripts/costs/total_maintenance_costs.py

@@ -57,14 +57,15 @@ class TotalMaintenanceCosts(CostBase):
     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.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 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:

scripts/costs/total_operational_costs.py

@@ -196,7 +196,7 @@ class TotalOperationalCosts(CostBase):
         if cooling is not None:
           hourly += cooling[i] / 3600
         if dhw is not None:
-          dhw += dhw[i] / 3600
+          hourly += dhw[i] / 3600
         hourly_fuel_consumption.append(hourly)
     else:
       heating = None

scripts/costs/total_operational_incomes.py

@@ -33,14 +33,12 @@ class TotalOperationalIncomes(CostBase):
       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)
-      # todo: check the adequate assignation of price. Pilar
-      price_export = archetype.income.electricity_export * cte.WATTS_HOUR_TO_JULES * 1000  # to account for unit change
+      price_export = archetype.income.electricity_export  # to account for unit change
       self._yearly_operational_incomes.loc[year, 'Incomes electricity'] = (
-          onsite_electricity_production * price_export * price_increase_electricity
+        (onsite_electricity_production / 3.6e6) * price_export * price_increase_electricity
       )
 
     self._yearly_operational_incomes.fillna(0, inplace=True)
-    return self._yearly_operational_incomes
+    return self._yearly_operational_incomes

scripts/district_heating_network/district_heating_network_creator.py

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

scripts/district_heating_network/geojson_graph_creator.py

@@ -0,0 +1,54 @@
+import json
+from shapely import LineString, Point
+import networkx as nx
+from pathlib import Path
+
+
+def networkx_to_geojson(graph: nx.Graph) -> Path:
+    """
+    Convert a NetworkX graph to GeoJSON format.
+
+    :param graph: A NetworkX graph.
+    :return: GeoJSON formatted dictionary.
+    """
+    features = []
+
+    for u, v, data in graph.edges(data=True):
+        line = LineString([u, v])
+        feature = {
+            "type": "Feature",
+            "geometry": {
+                "type": "LineString",
+                "coordinates": list(line.coords)
+            },
+            "properties": {
+                "weight": data.get("weight", 1.0)
+            }
+        }
+        features.append(feature)
+
+    for node, data in graph.nodes(data=True):
+        point = Point(node)
+        feature = {
+            "type": "Feature",
+            "geometry": {
+                "type": "Point",
+                "coordinates": list(point.coords)[0]
+            },
+            "properties": {
+                "type": data.get("type", "unknown"),
+                "id": data.get("id", "N/A")
+            }
+        }
+        features.append(feature)
+
+    geojson = {
+        "type": "FeatureCollection",
+        "features": features
+    }
+
+    output_geojson_file = Path('./out_files/network_graph.geojson').resolve()
+    with open(output_geojson_file, 'w') as file:
+        json.dump(geojson, file, indent=4)
+
+    return output_geojson_file

scripts/district_heating_network/road_processor.py

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

scripts/district_heating_network/simultinity_factor.py

@@ -0,0 +1,80 @@
+import pandas as pd
+import numpy as np
+
+class DemandShiftProcessor:
+  def __init__(self, city):
+    self.city = city
+
+  def random_shift(self, series):
+    shift_amount = np.random.randint(0, round(0.005 * len(series)))
+    return series.shift(shift_amount).fillna(series.shift(shift_amount - len(series)))
+
+  def process_demands(self):
+    heating_dfs = []
+    cooling_dfs = []
+
+    for building in self.city.buildings:
+      heating_df = self.convert_building_to_dataframe(building, 'heating')
+      cooling_df = self.convert_building_to_dataframe(building, 'cooling')
+      heating_df.set_index('Date/Time', inplace=True)
+      cooling_df.set_index('Date/Time', inplace=True)
+      shifted_heating_demands = heating_df.apply(self.random_shift, axis=0)
+      shifted_cooling_demands = cooling_df.apply(self.random_shift, axis=0)
+      self.update_building_demands(building, shifted_heating_demands, 'heating')
+      self.update_building_demands(building, shifted_cooling_demands, 'cooling')
+      heating_dfs.append(shifted_heating_demands)
+      cooling_dfs.append(shifted_cooling_demands)
+
+    combined_heating_df = pd.concat(heating_dfs, axis=1)
+    combined_cooling_df = pd.concat(cooling_dfs, axis=1)
+    self.calculate_and_set_simultaneity_factor(combined_heating_df, 'heating')
+    self.calculate_and_set_simultaneity_factor(combined_cooling_df, 'cooling')
+
+  def convert_building_to_dataframe(self, building, demand_type):
+    if demand_type == 'heating':
+      data = {
+        "Date/Time": self.generate_date_time_index(),
+        "Heating_Demand": building.heating_demand["hour"]
+      }
+    else:  # cooling
+      data = {
+        "Date/Time": self.generate_date_time_index(),
+        "Cooling_Demand": building.cooling_demand["hour"]
+      }
+    return pd.DataFrame(data)
+
+  def generate_date_time_index(self):
+    # Generate hourly date time index for a full year in 2024
+    date_range = pd.date_range(start="2013-01-01 00:00:00", end="2013-12-31 23:00:00", freq='H')
+    return date_range.strftime('%m/%d %H:%M:%S').tolist()
+
+  def update_building_demands(self, building, shifted_demands, demand_type):
+    if demand_type == 'heating':
+      shifted_series = shifted_demands["Heating_Demand"]
+      building.heating_demand = self.calculate_new_demands(shifted_series)
+    else:  # cooling
+      shifted_series = shifted_demands["Cooling_Demand"]
+      building.cooling_demand = self.calculate_new_demands(shifted_series)
+
+  def calculate_new_demands(self, shifted_series):
+    new_demand = {
+      "hour": shifted_series.tolist(),
+      "month": self.calculate_monthly_demand(shifted_series),
+      "year": [shifted_series.sum()]
+    }
+    return new_demand
+
+  def calculate_monthly_demand(self, series):
+    series.index = pd.to_datetime(series.index, format='%m/%d %H:%M:%S')
+    monthly_demand = series.resample('M').sum()
+    return monthly_demand.tolist()
+
+  def calculate_and_set_simultaneity_factor(self, combined_df, demand_type):
+    total_demand_original = combined_df.sum(axis=1)
+    peak_total_demand_original = total_demand_original.max()
+    individual_peak_demands = combined_df.max(axis=0)
+    sum_individual_peak_demands = individual_peak_demands.sum()
+    if demand_type == 'heating':
+      self.city.simultaneity_factor_heating = peak_total_demand_original / sum_individual_peak_demands
+    else:  # cooling
+      self.city.simultaneity_factor_cooling = peak_total_demand_original / sum_individual_peak_demands

scripts/energy_system_analysis_report.py

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

scripts/energy_system_retrofit_report.py

@@ -0,0 +1,596 @@
+import os
+import hub.helpers.constants as cte
+import matplotlib.pyplot as plt
+from matplotlib import cm
+from scripts.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 / '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']:
+        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] /
+                             (cte.WATTS_HOUR_TO_JULES * 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()

scripts/energy_system_retrofit_results.py

@@ -1,68 +1,176 @@
 import hub.helpers.constants as cte
 
 
-def system_results(buildings):
-  system_performance_summary = {}
-  fields = ["Energy System Archetype", "Heating Equipments", "Cooling Equipments", "DHW Equipments",
-            "Photovoltaic System Capacity", "Heating Fuel", "Yearly HVAC Energy Consumption (MWh)",
-            "DHW Energy Consumption (MWH)", "PV Yearly Production (kWh)", "LCC Analysis Duration (Years)",
-            "Energy System Capital Cost (CAD)", "Energy System Average Yearly Operational Cost (CAD)",
-            "Energy System Life Cycle Cost (CAD)"]
-  for building in buildings:
-    system_performance_summary[f'{building.name}'] = {}
-    for field in fields:
-      system_performance_summary[f'{building.name}'][field] = '-'
-
-  for building in buildings:
-    fuels = []
-    system_performance_summary[f'{building.name}']['Energy System Archetype'] = building.energy_systems_archetype_name
-    energy_systems = building.energy_systems
+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:
-      demand_types = energy_system.demand_types
-      for demand_type in demand_types:
-        if demand_type == cte.COOLING:
-          equipments = []
-          for generation_system in energy_system.generation_systems:
-            if generation_system.fuel_type == cte.ELECTRICITY:
-              equipments.append(generation_system.name or generation_system.system_type)
-            cooling_equipments = ", ".join(equipments)
-            system_performance_summary[f'{building.name}']['Cooling Equipments'] = cooling_equipments
-        elif demand_type == cte.HEATING:
-          equipments = []
-          for generation_system in energy_system.generation_systems:
-            if generation_system.nominal_heat_output is not None:
-              equipments.append(generation_system.name or generation_system.system_type)
-              fuels.append(generation_system.fuel_type)
-            heating_equipments = ", ".join(equipments)
-            system_performance_summary[f'{building.name}']['Heating Equipments'] = heating_equipments
-        elif demand_type == cte.DOMESTIC_HOT_WATER:
-          equipments = []
-          for generation_system in energy_system.generation_systems:
-            equipments.append(generation_system.name or generation_system.system_type)
-            dhw_equipments = ", ".join(equipments)
-            system_performance_summary[f'{building.name}']['DHW Equipments'] = dhw_equipments
-      for generation_system in energy_system.generation_systems:
-        if generation_system.system_type == cte.PHOTOVOLTAIC:
-          system_performance_summary[f'{building.name}'][
-            'Photovoltaic System Capacity'] = generation_system.nominal_electricity_output or str(0)
-    heating_fuels = ", ".join(fuels)
-    system_performance_summary[f'{building.name}']['Heating Fuel'] = heating_fuels
-    system_performance_summary[f'{building.name}']['Yearly HVAC Energy Consumption (MWh)'] = format(
-      (building.heating_consumption[cte.YEAR][0] + building.cooling_consumption[cte.YEAR][0]) / 3.6e9, '.2f')
-    system_performance_summary[f'{building.name}']['DHW Energy Consumption (MWH)'] = format(
-      building.domestic_hot_water_consumption[cte.YEAR][0] / 1e6, '.2f')
-  return system_performance_summary
+      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 new_system_results(buildings):
-  new_system_performance_summary = {}
-  fields = ["Energy System Archetype", "Heating Equipments", "Cooling Equipments", "DHW Equipments",
-            "Photovoltaic System Capacity", "Heating Fuel", "Yearly HVAC Energy Consumption (MWh)",
-            "DHW Energy Consumption (MWH)", "PV Yearly Production (kWh)", "LCC Analysis Duration (Years)",
-            "Energy System Capital Cost (CAD)", "Energy System Average Yearly Operational Cost (CAD)",
-            "Energy System Life Cycle Cost (CAD)"]
-  for building in buildings:
-    new_system_performance_summary[f'{building.name}'] = {}
-    for field in fields:
-      new_system_performance_summary[f'{building.name}'][field] = '-'
-  return new_system_performance_summary
+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

scripts/energy_system_sizing_and_simulation_factory.py

@@ -8,6 +8,7 @@ Project Coder Saeed Ranjbar saeed.ranjbar@mail.concordia.ca
 from scripts.system_simulation_models.archetype13 import Archetype13
 from scripts.system_simulation_models.archetype13_stratified_tes import Archetype13Stratified
 from scripts.system_simulation_models.archetype1 import Archetype1
+from scripts.system_simulation_models.archetypes14_15 import Archetype14_15
 
 
 class EnergySystemsSimulationFactory:
@@ -36,6 +37,15 @@ class EnergySystemsSimulationFactory:
     self._building.level_of_detail.energy_systems = 2
     self._building.level_of_detail.energy_systems = 2
 
+  def _archetype14_15(self):
+    """
+    Enrich the city by using the sizing and simulation model developed for archetype14 and archetype15 of
+    montreal_future_systems
+    """
+    Archetype14_15(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

scripts/ep_run_enrich.py

@@ -9,10 +9,10 @@ from hub.imports.results_factory import ResultFactory
 sys.path.append('./')
 
 
-def energy_plus_workflow(city):
+def energy_plus_workflow(city, output_path):
   try:
     # city = city
-    out_path = (Path(__file__).parent.parent / 'out_files')
+    out_path = output_path
     files = glob.glob(f'{out_path}/*')
 
     # for file in files:

scripts/geojson_creator.py

@@ -4,13 +4,16 @@ from shapely import Point
 from pathlib import Path
 
 
-def process_geojson(x, y, diff):
+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()
-  output_file = Path('./input_files/output_buildings.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:

scripts/pv_feasibility.py

@@ -0,0 +1,37 @@
+from pathlib import Path
+import subprocess
+from hub.imports.geometry_factory import GeometryFactory
+from scripts.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 / 'input_files')
+  output_path = (Path(__file__).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
+
+
+
+

scripts/pv_sizing_and_simulation.py

@@ -49,18 +49,11 @@ class PVSizingSimulation(RadiationTilted):
     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 for x in radiation]
+    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])]
-
-
-
-
-
-
-
-
-
+    self.building.onsite_electrical_production[cte.YEAR] = [sum(self.building.onsite_electrical_production[cte.MONTH])]

scripts/random_assignation.py

@@ -15,8 +15,8 @@ from hub.city_model_structure.building import Building
 energy_systems_format = 'montreal_custom'
 
 # parameters:
-residential_systems_percentage = {'system 1 gas': 100,
-                                  'system 1 electricity': 0,
+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,
@@ -25,11 +25,11 @@ residential_systems_percentage = {'system 1 gas': 100,
                                   'system 5 electricity': 0,
                                   'system 6 gas': 0,
                                   'system 6 electricity': 0,
-                                  'system 8 gas': 0,
-                                  'system 8 electricity': 0}
+                                  'system 8 gas': 15,
+                                  'system 8 electricity': 35}
 
-residential_new_systems_percentage = {'PV+ASHP+GasBoiler+TES': 0,
-                                      'PV+4Pipe+DHW': 100,
+residential_new_systems_percentage = {'PV+ASHP+GasBoiler+TES': 100,
+                                      'PV+4Pipe+DHW': 0,
                                       'PV+ASHP+ElectricBoiler+TES': 0,
                                       'PV+GSHP+GasBoiler+TES': 0,
                                       'PV+GSHP+ElectricBoiler+TES': 0,

scripts/report_creation.py

@@ -1,73 +1,119 @@
 import subprocess
 import datetime
+import os
+from pathlib import Path
+
 
 class LatexReport:
-  def __init__(self, file_name):
-      self.file_name = file_name
-      self.content = []
-      self.content.append(r'\documentclass{article}')
-      self.content.append(r'\usepackage[margin=2.5cm]{geometry}')  # Adjust page margins
-      self.content.append(r'\usepackage{graphicx}')
-      self.content.append(r'\usepackage{tabularx}')
-      self.content.append(r'\begin{document}')
-      # Get current date and time
-      current_datetime = datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S")
-      self.content.append(r'\title{Energy System Analysis Report - ' + current_datetime + r'}')
-      self.content.append(r'\author{Next-Generation Cities Institute}')
-      self.content.append(r'\date{}')  # Remove the date field, as it's included in the title now
-      self.content.append(r'\maketitle')
+  def __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'}')
+    self.content.append(r'\section{' + section_title + r'}')
 
   def add_subsection(self, subsection_title):
-      self.content.append(r'\subsection{' + subsection_title + r'}')
+    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)
+    self.content.append(text)
 
-  def add_table(self, table_data, caption=None, first_column_width=None):
+  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  # Default total width
+    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  # Adjust total width for the first column
+      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')
 
-    self.content.append(r'\begin{tabularx}{\textwidth}{|p{' + first_column_width_str + r'}|' + '|'.join(['X'] * (
-          num_columns - 1)) + '|}' if first_column_width is not None else r'\begin{tabularx}{\textwidth}{|' + '|'.join(
-      ['X'] * num_columns) + '|}')
+    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')
-    for row in table_data:
-      self.content.append(' & '.join(row) + r' \\')
+
+    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):
+  def add_image(self, image_path, caption=None, placement='ht'):
     if caption:
-      self.content.append(r'\begin{figure}[htbp]')
+      self.content.append(r'\begin{figure}[' + placement + r']')
       self.content.append(r'\centering')
-      self.content.append(r'\includegraphics[width=0.8\textwidth]{' + image_path + r'}')
+      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}[htbp]')
+      self.content.append(r'\begin{figure}[' + placement + r']')
       self.content.append(r'\centering')
-      self.content.append(r'\includegraphics[width=0.8\textwidth]{' + image_path + r'}')
+      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}')  # Add this line to close the document
-    with open(self.file_name, 'w') as f:
+    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', self.file_name])
+    subprocess.run(['pdflatex', '-output-directory', str(self.output_path), str(self.file_path)])
+

scripts/system_simulation.py

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

scripts/system_simulation_models/archetype1.py

@@ -19,6 +19,8 @@ class Archetype1:
     self._hourly_dhw_demand = building.domestic_hot_water_heat_demand[cte.HOUR]
     self._output_path = output_path
     self._t_out = [0] + building.external_temperature[cte.HOUR]
+    self.results = {}
+    self.dt = 900
 
   def hvac_sizing(self):
     storage_factor = 3
@@ -29,13 +31,26 @@ class Archetype1:
     heat_pump.nominal_cooling_output = round(self._cooling_peak_load / 3600)
     boiler.nominal_heat_output = round(0.5 * self._heating_peak_load / 3600)
     thermal_storage.volume = round(
-      (self._heating_peak_load * storage_factor) / (cte.WATER_HEAT_CAPACITY * cte.WATER_DENSITY * 30))
+      (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))
+    return dhw_hp, dhw_tes
+
   def heating_system_simulation(self):
     hp, boiler, tes = self.hvac_sizing()
     cop_curve_coefficients = [float(coefficient) for coefficient in hp.heat_efficiency_curve.coefficients]
-    demand = self._hourly_heating_demand
+    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
     # Heating System Simulation
     variable_names = ["t_sup_hp", "t_tank", "t_ret", "m_ch", "m_dis", "q_hp", "q_boiler", "hp_cop",
@@ -61,7 +76,7 @@ class Archetype1:
     for i in range(len(demand) - 1):
       t_tank[i + 1] = (t_tank[i] +
                        ((m_ch[i] * (t_sup_hp[i] - t_tank[i])) +
-                        (ua * (self._t_out[i] - t_tank[i] + 5)) / cte.WATER_HEAT_CAPACITY -
+                        (ua * (t_out[i] - t_tank[i])) / cte.WATER_HEAT_CAPACITY -
                         m_dis[i] * (t_tank[i] - t_ret[i])) * (dt / (cte.WATER_DENSITY * v)))
       if t_tank[i + 1] < 40:
         q_hp[i + 1] = hp_heating_cap
@@ -77,16 +92,16 @@ class Archetype1:
         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 * self._t_out[i + 1] + 32
+      t_tank_fahrenheit = 1.8 * t_tank[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 +
+        hp_cop[i + 1] = (1 / (cop_curve_coefficients[0] +
+                         cop_curve_coefficients[1] * t_tank_fahrenheit +
+                         cop_curve_coefficients[2] * t_tank_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] / hp_efficiency
+                         cop_curve_coefficients[5] * t_tank_fahrenheit * t_out_fahrenheit)) * hp_efficiency
+        hp_electricity[i + 1] = q_hp[i + 1] / hp_cop[i + 1]
       else:
         hp_cop[i + 1] = 0
         hp_electricity[i + 1] = 0
@@ -110,101 +125,253 @@ class Archetype1:
       boiler_gas[i + 1] = (q_boiler[i + 1] * dt) / cte.NATURAL_GAS_LHV
       boiler_consumption[i + 1] = q_boiler[i + 1] / boiler_efficiency
       heating_consumption[i + 1] = boiler_consumption[i + 1] + hp_electricity[i + 1]
-    data = list(zip(self._t_out, q_hp, hp_electricity, hp_cop, m_ch, t_sup_hp, t_tank, m_dis, q_boiler,
-                    boiler_gas, t_sup_boiler, t_ret, demand))
-    file_name = f'heating_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(['Outside_air_temperature', 'HP_heat_output(W)', 'HP_heating_electricity_consumption(W)',
-                            'HP_COP', 'TES_charge_flow_rate(kg/s)', 'TES_supply_temperature(C)', 'TES_temperature(C)',
-                            'TES_discharge_flow_rate(kg/s)', 'Boiler_heat_output(W)', 'Boiler_gas_consumption(m3)',
-                            'Boiler_supply_temperature(C)', 'Heating_loop_return_temperature(C)', 'Heating_demand (W)'])
-      # Write data
-      output_file.writerows(data)
-    return heating_consumption, hp_electricity, boiler_consumption, t_sup_hp, t_sup_boiler
+    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_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_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)
-    demand = self._hourly_cooling_demand
+    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"]
+    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_eer) = [variables[name] 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] = 13
-    dt = 3600
-    for i in range(len(demand) - 1):
-      if demand[i] > 0:
-        m[i] = self._cooling_peak_load / (cte.WATER_HEAT_CAPACITY * 5 * dt)
-        if t_ret[i] > 13:
-          q_hp[i] = hp.nominal_cooling_output
-          t_sup_hp[i] = t_ret[i] - q_hp[i] / (m[i] * cte.WATER_HEAT_CAPACITY)
+
+    for i in range(1, len(demand)):
+      if demand[i] > 0.15 * self._cooling_peak_load:
+        m[i] = hp.nominal_cooling_output / (cte.WATER_HEAT_CAPACITY * 5)
+        if t_ret[i - 1] >= 13:
+          if demand[i] < 0.25 * self._cooling_peak_load:
+            q_hp[i] = 0.25 * hp.nominal_cooling_output
+          elif demand[i] < 0.5 * self._cooling_peak_load:
+            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]
-        t_ret[i + 1] = t_sup_hp[i] + demand[i] / (m[i] * cte.WATER_HEAT_CAPACITY)
+          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]
-        t_ret[i + 1] = t_ret[i]
+        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 * self._t_out[i] + 32
+      t_out_fahrenheit = 1.8 * t_out[i] + 32
       if q_hp[i] > 0:
-        hp_eer[i] = 1/((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_cop[i] = (1 / (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)) * cooling_efficiency / 3.41
         hp_electricity[i] = q_hp[i] / cooling_efficiency
       else:
-        hp_eer[i] = 0
+        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_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_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 hp_hourly
+
+  def dhw_system_simulation(self):
+    hp, tes = self.dhw_sizing()
+    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] < 62:
+        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))
+      if t_tank[i] < 60:
+        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_out[i] +
+                     cop_curve_coefficients[2] * t_out[i] ** 2 +
+                     cop_curve_coefficients[3] * t_tank[i] +
+                     cop_curve_coefficients[4] * t_tank[i] ** 2 +
+                     cop_curve_coefficients[5] * t_tank[i] * t_out[i]) * float(hp.heat_efficiency)
+        hp_electricity[i] = q_hp[i] / hp_cop[i]
+      else:
+        hp_cop[i] = 0
         hp_electricity[i] = 0
 
-    data = list(zip(q_hp, hp_electricity, hp_eer, m, t_sup_hp, t_ret, demand, self._t_out))
-    file_name = f'cooling_system_simulation_results_{self._name}.csv'
+      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(['HP_cooling_output (W)', 'HP_cooling_electricity_consumption (W)', 'HP_EER',
-                            'Cooling_loop_flow_rate(kg/s)', 'Cooling_supply_temperature(C)',
-                            'Cooling_loop_return_temperature(C)', 'Cooling_demand (W)', 'Outside_Air_Temperature'])
+      output_file.writerow(self.results.keys())
       # Write data
-      output_file.writerows(data)
-      return hp_electricity, t_sup_hp
+      output_file.writerows(zip(*self.results.values()))
 
-
-  def enrich_buildings(self):
-    (self._building.heating_consumption[cte.HOUR], hp_electricity_heating,
-     boiler_consumption, t_sup_hp, t_sup_boiler) = self.heating_system_simulation()
-    hp_electricity_cooling, t_sup_hp_cooling = self.cooling_system_simulation()
-    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._hvac_system.generation_systems[0].energy_consumption[cte.HEATING] = {}
-    self._hvac_system.generation_systems[1].energy_consumption[cte.HEATING] = {}
-    self._hvac_system.generation_systems[0].energy_consumption[cte.HEATING][cte.HOUR] = hp_electricity_heating
-    self._hvac_system.generation_systems[0].energy_consumption[cte.HEATING][cte.MONTH] = MonthlyValues.get_total_month(
-      self._hvac_system.generation_systems[0].energy_consumption[cte.HEATING][cte.HOUR])
-    self._hvac_system.generation_systems[0].energy_consumption[cte.HEATING][cte.YEAR] = \
-        [sum(self._hvac_system.generation_systems[0].energy_consumption[cte.HEATING][cte.MONTH])]
-    self._hvac_system.generation_systems[1].energy_consumption[cte.HEATING][cte.HOUR] = boiler_consumption
-    self._hvac_system.generation_systems[1].energy_consumption[cte.HEATING][cte.MONTH] = MonthlyValues.get_total_month(
-      self._hvac_system.generation_systems[1].energy_consumption[cte.HEATING][cte.HOUR])
-    self._hvac_system.generation_systems[1].energy_consumption[cte.HEATING][cte.YEAR] = \
-        [sum(self._hvac_system.generation_systems[1].energy_consumption[cte.HEATING][cte.MONTH])]
-    self._hvac_system.generation_systems[0].heat_supply_temperature = t_sup_hp
-    self._hvac_system.generation_systems[1].heat_supply_temperature = t_sup_boiler
-    self._hvac_system.generation_systems[0].energy_consumption[cte.COOLING] = {}
-    self._hvac_system.generation_systems[0].energy_consumption[cte.COOLING][cte.HOUR] = hp_electricity_cooling
-    self._hvac_system.generation_systems[0].energy_consumption[cte.COOLING][cte.MONTH] = MonthlyValues.get_total_month(
-      self._hvac_system.generation_systems[0].energy_consumption[cte.COOLING][cte.HOUR])
-    self._hvac_system.generation_systems[0].energy_consumption[cte.COOLING][cte.YEAR] = \
-        [sum(self._hvac_system.generation_systems[0].energy_consumption[cte.COOLING][cte.MONTH])]
-    self._hvac_system.generation_systems[0].cooling_supply_temperature = t_sup_hp_cooling
-    self._dhw_system.generation_systems[0].energy_consumption[cte.DOMESTIC_HOT_WATER] = (
-      self._building.domestic_hot_water_consumption)

scripts/system_simulation_models/archetype13.py

@@ -17,9 +17,10 @@ class Archetype13:
     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 = building.domestic_hot_water_heat_demand[cte.HOUR]
+    self._hourly_heating_demand = [demand / cte.WATTS_HOUR_TO_JULES for demand in building.heating_demand[cte.HOUR]]
+    self._hourly_cooling_demand = [demand / cte.WATTS_HOUR_TO_JULES 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 = {}
@@ -30,11 +31,12 @@ class Archetype13:
     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(0.5 * self._heating_peak_load / 3600)
-    heat_pump.nominal_cooling_output = round(self._cooling_peak_load / 3600)
-    boiler.nominal_heat_output = round(0.5 * self._heating_peak_load / 3600)
+    heat_pump.nominal_heat_output = round(0.5 * self._heating_peak_load)
+    heat_pump.nominal_cooling_output = round(self._cooling_peak_load)
+    boiler.nominal_heat_output = round(0.5 * self._heating_peak_load)
     thermal_storage.volume = round(
-      (self._heating_peak_load * storage_factor) / (cte.WATER_HEAT_CAPACITY * cte.WATER_DENSITY * 25))
+      (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):
@@ -64,6 +66,7 @@ class Archetype13:
         [variables[name] for name in variable_names]
     t_tank[0] = 55
     hp_heating_cap = hp.nominal_heat_output
+    hp_efficiency = float(hp.heat_efficiency)
     boiler_heating_cap = boiler.nominal_heat_output
     hp_delta_t = 5
     boiler_efficiency = float(boiler.heat_efficiency)
@@ -96,15 +99,15 @@ class Archetype13:
         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_tank_fahrenheit = 1.8 * t_tank[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 +
+        hp_cop[i + 1] = (1 / (cop_curve_coefficients[0] +
+                         cop_curve_coefficients[1] * t_tank_fahrenheit +
+                         cop_curve_coefficients[2] * t_tank_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)
+                         cop_curve_coefficients[5] * t_tank_fahrenheit * t_out_fahrenheit)) * hp_efficiency
         hp_electricity[i + 1] = q_hp[i + 1] / hp_cop[i + 1]
       else:
         hp_cop[i + 1] = 0
@@ -123,11 +126,11 @@ class Archetype13:
         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:
+        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 * cte.HOUR_TO_SECONDS)
+        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)
     tes.temperature = []
     hp_electricity_j = [(x * cte.WATTS_HOUR_TO_JULES) / number_of_ts for x in hp_electricity]
@@ -181,19 +184,19 @@ class Archetype13:
     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"]
+    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_eer) = [variables[name] 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] = 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 demand[i] > 0.15 * self._cooling_peak_load:
+        m[i] = hp.nominal_cooling_output / (cte.WATER_HEAT_CAPACITY * 5)
         if t_ret[i - 1] >= 13:
-          if demand[i] < 0.25 * self._cooling_peak_load / cte.HOUR_TO_SECONDS:
+          if demand[i] < 0.25 * self._cooling_peak_load:
             q_hp[i] = 0.25 * hp.nominal_cooling_output
-          elif demand[i] < 0.5 * self._cooling_peak_load / cte.HOUR_TO_SECONDS:
+          elif demand[i] < 0.5 * self._cooling_peak_load:
             q_hp[i] = 0.5 * hp.nominal_cooling_output
           else:
             q_hp[i] = hp.nominal_cooling_output
@@ -208,20 +211,20 @@ class Archetype13:
       else:
         m[i] = 0
         q_hp[i] = 0
-        t_sup_hp[i] = t_ret[i -1]
+        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] +
+        hp_cop[i] = (1 / (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] / hp_eer[i]
+                     eer_curve_coefficients[5] * t_sup_hp_fahrenheit * t_out_fahrenheit)) * cooling_efficiency / 3.41
+        hp_electricity[i] = q_hp[i] / cooling_efficiency
       else:
-        hp_eer[i] = 0
+        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 = []
@@ -240,7 +243,7 @@ class Archetype13:
     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 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
@@ -274,7 +277,7 @@ class Archetype13:
     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:
+      if t_tank[i] < 62:
         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:
@@ -282,14 +285,8 @@ class Archetype13:
         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)
+      if t_tank[i] < 60:
+        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:
@@ -302,11 +299,11 @@ class Archetype13:
       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)
+                     cop_curve_coefficients[1] * t_out[i] +
+                     cop_curve_coefficients[2] * t_out[i] ** 2 +
+                     cop_curve_coefficients[3] * t_tank[i] +
+                     cop_curve_coefficients[4] * t_tank[i] ** 2 +
+                     cop_curve_coefficients[5] * t_tank[i] * t_out[i]) * float(hp.heat_efficiency)
         hp_electricity[i] = q_hp[i] / hp_cop[i]
       else:
         hp_cop[i] = 0
@@ -375,8 +372,8 @@ class Archetype13:
     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]))
+    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)

scripts/system_simulation_models/archetype13_stratified_tes.py

@@ -24,6 +24,7 @@ class Archetype13Stratified:
     self._output_path = output_path
     self._t_out = building.external_temperature[cte.HOUR]
     self.results = {}
+    self.dt = 300
 
   def hvac_sizing(self):
     storage_factor = 3
@@ -49,6 +50,7 @@ class Archetype13Stratified:
 
   def heating_system_simulation_stratified(self):
     hp, boiler, tes = self.hvac_sizing()
+    hp_efficiency = float(hp.heat_efficiency)
     cop_curve_coefficients = [float(coefficient) for coefficient in hp.heat_efficiency_curve.coefficients]
     demand = [0] + [x for x in self._hourly_heating_demand for _ in range(12)]
     hp.source_temperature = self._t_out
@@ -114,15 +116,15 @@ class Archetype13Stratified:
         t_sup_hp[i + 1] = (q_hp[i + 1] / (m_ch[i + 1] * cte.WATER_HEAT_CAPACITY)) + t4[i + 1]
       else:
         q_hp[i + 1], m_ch[i + 1], t_sup_hp[i + 1] = 0, 0, t4[i + 1]
-      t_sup_hp_fahrenheit = 1.8 * t_sup_hp[i + 1] + 32
+      t_tank_fahrenheit = 1.8 * t4[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 +
+        hp_cop[i + 1] = (1 / (cop_curve_coefficients[0] +
+                         cop_curve_coefficients[1] * t_tank_fahrenheit +
+                         cop_curve_coefficients[2] * t_tank_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)
+                         cop_curve_coefficients[5] * t_tank_fahrenheit * t_out_fahrenheit)) * hp_efficiency
         hp_electricity[i + 1] = q_hp[i + 1] / hp_cop[i + 1]
       else:
         hp_cop[i + 1] = 0
@@ -211,66 +213,81 @@ class Archetype13Stratified:
     hp = self.hvac_sizing()[0]
     eer_curve_coefficients = [float(coefficient) for coefficient in hp.cooling_efficiency_curve.coefficients]
     cooling_efficiency = float(hp.cooling_efficiency)
-    demand = self._hourly_cooling_demand
+    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"]
+    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_eer) = [variables[name] 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] = 13
-    dt = 3600
-    for i in range(len(demand) - 1):
-      if demand[i] > 0:
-        m[i] = self._cooling_peak_load / (cte.WATER_HEAT_CAPACITY * 5 * dt)
-        if t_ret[i] > 13:
-          if demand[i] < 0.25 * self._cooling_peak_load / dt:
+
+    for i in range(1, len(demand)):
+      if demand[i] > 0.15 * self._cooling_peak_load:
+        m[i] = hp.nominal_cooling_output / (cte.WATER_HEAT_CAPACITY * 5)
+        if t_ret[i - 1] >= 13:
+          if demand[i] < 0.25 * self._cooling_peak_load:
             q_hp[i] = 0.25 * hp.nominal_cooling_output
-          elif demand[i] < 0.5 * self._cooling_peak_load / dt:
+          elif demand[i] < 0.5 * self._cooling_peak_load:
             q_hp[i] = 0.5 * hp.nominal_cooling_output
           else:
             q_hp[i] = hp.nominal_cooling_output
-          t_sup_hp[i] = t_ret[i] - q_hp[i] / (m[i] * cte.WATER_HEAT_CAPACITY)
+          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]
-        t_ret[i + 1] = t_sup_hp[i] + demand[i] / (m[i] * cte.WATER_HEAT_CAPACITY)
+          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]
-        t_ret[i + 1] = t_ret[i]
+        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 * self._t_out[i] + 32
+      t_out_fahrenheit = 1.8 * t_out[i] + 32
       if q_hp[i] > 0:
-        hp_eer[i] = (eer_curve_coefficients[0] +
+        hp_cop[i] = (1 / (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)
+                     eer_curve_coefficients[5] * t_sup_hp_fahrenheit * t_out_fahrenheit)) * cooling_efficiency / 3.41
         hp_electricity[i] = q_hp[i] / cooling_efficiency
       else:
-        hp_eer[i] = 0
+        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_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_electricity
+    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_electricity
+    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_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 hp_hourly
 
   def dhw_system_simulation(self):
     hp, tes = self.dhw_sizing()
     cop_curve_coefficients = [float(coefficient) for coefficient in hp.heat_efficiency_curve.coefficients]
-    demand = self._hourly_dhw_demand
+    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)
@@ -279,7 +296,7 @@ class Archetype13Stratified:
         [variables[name] for name in variable_names]
     t_tank[0] = 70
     v_dhw[0] = tes.volume
-    dt = 3600
+
     hp_heating_cap = hp.nominal_heat_output
     hp_delta_t = 8
     v, h = float(tes.volume), float(tes.height)
@@ -292,26 +309,18 @@ class Archetype13Stratified:
     ua = u_tot * (2 * a_top + a_side)
     freshwater_temperature = 18
     for i in range(len(demand) - 1):
-      delta_t_demand = demand[i] * (dt / (cte.WATER_DENSITY * cte.WATER_HEAT_CAPACITY * v))
-      if t_tank[i] < 65:
+      delta_t_demand = demand[i] * (self.dt / (cte.WATER_DENSITY * cte.WATER_HEAT_CAPACITY * v))
+      if t_tank[i] < 62:
         q_hp[i] = hp_heating_cap
-      delta_t_hp = q_hp[i] * (dt / (cte.WATER_DENSITY * cte.WATER_HEAT_CAPACITY * v))
+      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) * (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) / dt
-            if power <= float(tes.heating_coil_capacity):
-              q_coil[i] = power
-            else:
-              q_coil[i] = float(tes.heating_coil_capacity)
-          elif t_tank[i] < 65:
-            q_coil[i] = float(tes.heating_coil_capacity)
-      delta_t_coil = q_coil[i] * (dt / (cte.WATER_DENSITY * cte.WATER_HEAT_CAPACITY * v))
+      delta_t_freshwater = m_refill[i] * (t_tank[i] - freshwater_temperature) * (self.dt / (v * cte.WATER_DENSITY))
+      if t_tank[i] < 60:
+        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)
@@ -320,69 +329,84 @@ class Archetype13Stratified:
         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 * self._t_out[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] / 3.5
+                     cop_curve_coefficients[1] * t_out[i] +
+                     cop_curve_coefficients[2] * t_out[i] ** 2 +
+                     cop_curve_coefficients[3] * t_tank[i] +
+                     cop_curve_coefficients[4] * t_tank[i] ** 2 +
+                     cop_curve_coefficients[5] * t_tank[i] * t_out[i]) * float(hp.heat_efficiency)
+        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)
+    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_electricity
+    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] = q_coil
+    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'] = self._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_electricity, q_coil
+    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_stratified()
-    # hp_cooling = self.cooling_system_simulation()
-    # hp_dhw, heating_coil = self.dhw_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]))
+    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)

scripts/system_simulation_models/archetypes14_15.py

@@ -0,0 +1,398 @@
+import math
+import hub.helpers.constants as cte
+import csv
+from hub.helpers.monthly_values import MonthlyValues
+
+
+class Archetype14_15:
+  def __init__(self, building, output_path):
+    self._building = building
+    self._name = building.name
+    if 'PV' in building.energy_systems_archetype_name:
+      i = 1
+      self._pv_system = building.energy_systems[0]
+    else:
+      i = 0
+    self._dhw_system = building.energy_systems[i]
+    self._heating_system = building.energy_systems[i + 1]
+    self._cooling_system = building.energy_systems[i + 2]
+    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.WATTS_HOUR_TO_JULES for demand in building.heating_demand[cte.HOUR]]
+    self._hourly_cooling_demand = [demand / cte.WATTS_HOUR_TO_JULES 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 heating_system_sizing(self):
+    storage_factor = 3
+    heat_pump = self._heating_system.generation_systems[1]
+    heat_pump.source_temperature = self._t_out
+    boiler = self._heating_system.generation_systems[0]
+    thermal_storage = boiler.energy_storage_systems[0]
+    heat_pump.nominal_heat_output = round(0.5 * self._heating_peak_load)
+    boiler.nominal_heat_output = round(0.5 * self._heating_peak_load)
+    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 cooling_system_sizing(self):
+    heat_pump = self._cooling_system.generation_systems[0]
+    heat_pump.nominal_cooling_output = heat_pump.nominal_cooling_output = round(self._cooling_peak_load)
+    heat_pump.source_temperature = self._t_out
+    return heat_pump
+
+
+  def dhw_system_sizing(self):
+    storage_factor = 3
+    dhw_hp = self._dhw_system.generation_systems[0]
+    dhw_hp.nominal_heat_output = round(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 * cte.WATTS_HOUR_TO_JULES) /
+      (cte.WATER_HEAT_CAPACITY * cte.WATER_DENSITY * 10))
+    return dhw_hp, dhw_tes
+
+  def heating_system_simulation(self):
+    hp, boiler, tes = self.heating_system_sizing()
+    hp_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 = 5
+    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_tank_fahrenheit = 1.8 * t_tank[i + 1] + 32
+      t_out_fahrenheit = 1.8 * t_out[i + 1] + 32
+      if q_hp[i + 1] > 0:
+        hp_cop[i + 1] = (1 / (cop_curve_coefficients[0] +
+                         cop_curve_coefficients[1] * t_tank_fahrenheit +
+                         cop_curve_coefficients[2] * t_tank_fahrenheit ** 2 +
+                         cop_curve_coefficients[3] * t_out_fahrenheit +
+                         cop_curve_coefficients[4] * t_out_fahrenheit ** 2 +
+                         cop_curve_coefficients[5] * t_tank_fahrenheit * t_out_fahrenheit)) * hp_efficiency
+        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_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.cooling_system_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_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] = 13
+
+    for i in range(1, len(demand)):
+      if demand[i] > 0.15 * self._cooling_peak_load:
+        m[i] = hp.nominal_cooling_output / (cte.WATER_HEAT_CAPACITY * 5)
+        if t_ret[i - 1] >= 13:
+          if demand[i] < 0.25 * self._cooling_peak_load:
+            q_hp[i] = 0.25 * hp.nominal_cooling_output
+          elif demand[i] < 0.5 * self._cooling_peak_load:
+            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_cop[i] = (1 / (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)) * cooling_efficiency / 3.41
+        hp_electricity[i] = q_hp[i] / cooling_efficiency
+      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_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_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 hp_hourly
+
+  def dhw_system_simulation(self):
+    hp, tes = self.dhw_system_sizing()
+    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] < 62:
+        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))
+      if t_tank[i] < 60:
+        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_out[i] +
+                     cop_curve_coefficients[2] * t_out[i] ** 2 +
+                     cop_curve_coefficients[3] * t_tank[i] +
+                     cop_curve_coefficients[4] * t_tank[i] ** 2 +
+                     cop_curve_coefficients[5] * t_tank[i] * t_out[i]) * float(hp.heat_efficiency)
+        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)
+    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()))

simulation_result_test.py

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

tests/test_systems_catalog.py

@@ -38,10 +38,10 @@ class TestSystemsCatalog(TestCase):
     catalog = EnergySystemsCatalogFactory('montreal_future').catalog
 
     catalog_categories = catalog.names()
-    archetypes = catalog.names('archetypes')
-    self.assertEqual(13, len(archetypes['archetypes']))
+    archetypes = catalog.names()
+    self.assertEqual(15, len(archetypes['archetypes']))
     systems = catalog.names('systems')
-    self.assertEqual(10, len(systems['systems']))
+    self.assertEqual(12, len(systems['systems']))
     generation_equipments = catalog.names('generation_equipments')
     self.assertEqual(27, len(generation_equipments['generation_equipments']))
     with self.assertRaises(ValueError):