20240808 final model, archetype 13 is changed by commenting the lines of the curve coefficients and to include as input t_out = 15; main processes 4 scenarios

clusters.py

@@ -0,0 +1,124 @@
+
+import numpy as np
+import pandas as pd
+from sklearn.preprocessing import StandardScaler
+from sklearn.decomposition import PCA
+from sklearn.cluster import KMeans
+import matplotlib.pyplot as plt
+
+
+def extract_building_data(scenario, demand_types):
+    """
+    Extracts energy demand data for each building from the scenario dictionary.
+
+    Args:
+        scenario (dict): Scenario dictionary containing building data.
+        demand_types (list): List of demand types to extract.
+
+    Returns:
+        pd.DataFrame: DataFrame with buildings as rows and demand types as columns.
+    """
+    data = []
+    building_ids = []
+
+    for building_id, demand_data in scenario.items():
+        building_row = []
+        for demand_type in demand_types:
+            building_row.extend(demand_data[demand_type])
+        data.append(building_row)
+        building_ids.append(building_id)
+
+    return pd.DataFrame(data, index=building_ids,
+                        columns=[f"{demand_type}_{i}" for demand_type in demand_types for i in range(8760)])
+
+# Plot the data
+print('test')
+def cluster_buildings(scenario, demand_types, n_clusters=4,n_components=2):
+    """
+    Clusters buildings based on their energy demand.
+
+    Args:
+        scenario (dict): Scenario dictionary containing building data.
+        demand_types (list): List of demand types to use for clustering.
+        n_clusters (int): Number of clusters to form.
+
+    Returns:
+        pd.DataFrame: DataFrame with building IDs and their corresponding cluster labels.
+    """
+    # Extract building data
+    building_data = extract_building_data(scenario, demand_types)
+
+    # Standardize the data
+    scaler = StandardScaler()
+    building_data_scaled = scaler.fit_transform(building_data)
+
+    # Optional: Dimensionality reduction with PCA
+    pca = PCA(n_components)
+    building_data_pca = pca.fit_transform(building_data_scaled)
+
+    # Clustering with K-Means
+    kmeans = KMeans(n_clusters=n_clusters, random_state=42)
+    kmeans.fit(building_data_pca)
+
+    # Assign cluster labels to buildings
+    clusters = kmeans.labels_
+
+    # Create a DataFrame with building IDs and cluster labels
+    clustered_buildings = pd.DataFrame({
+        'building_id': building_data.index,
+        'cluster': clusters
+    })
+
+    return clustered_buildings
+
+
+# Example usage
+scenario_data = scenario['efficient with PV+4Pipe+DHW']  # Use 'baseline' or 'efficient with PV'
+demand_types = [
+    'heating_consumption_kWh'
+]
+
+# Cluster buildings
+clustered_buildings = cluster_buildings(scenario_data, demand_types, n_clusters=4,n_components=20)
+
+# Visualize clusters
+plt.figure(figsize=(10, 6))
+plt.scatter(clustered_buildings['building_id'], clustered_buildings['cluster'], c=clustered_buildings['cluster'],
+            cmap='viridis')
+plt.xlabel('Building ID')
+plt.ylabel('Cluster')
+plt.title('Building Clusters Based on Energy Demand')
+plt.colorbar(label='Cluster')
+# Save the plot
+plt.savefig(os.path.join(output_path, 'clusters.png'))
+plt.close()
+
+# Plot the data
+print('test')
+
+
+
+output_path_clusters = output_path / 'clustered_buildings_4.xlsx'
+clustered_buildings.to_excel(output_path_clusters, index=False)
+
+#clusters made in QGIS
+clusters_path=output_path/ "clusters" / 'updated_buildings_with_clusters.geojson'
+
+with open(clusters_path, 'r') as f:
+    clusters = json.load(f)
+
+def extract_cluster(key, clusters):
+
+    for idx, feature in enumerate(clusters['features']):
+        if str(feature['properties']['id']) == str(key):
+            cluster=feature['properties']['Cluster_ID']
+
+    return cluster
+cluster_dic={}
+scenario_data = scenario['efficient with PV+4Pipe+DHW']
+for key, data in scenario_data.items():
+    cluster_dic[key]=extract_cluster(key, clusters)
+
+#
+# cluster_df = pd.DataFrame.from_dict(cluster_dic, orient='index')
+# cluster_df.to_csv(output_path/'clusters.csv')

comparing.py

@@ -52,14 +52,6 @@ percentage_data = {
     1690: {"type1_%": 0.94, "type2_%": 0.06, "type3_%": 0, "roof_area": 2189.732519, "total_floor_area": 28466.52275},
     1691: {"type1_%": 0.75, "type2_%": 0.25, "type3_%": 0, "roof_area": 3159.077523, "total_floor_area": 31590.77523},
 }
-# Define the demand types to be plotted
-demand_types = [
-    'heating_demand_kWh',
-    'cooling_demand_kWh',
-    'domestic_hot_water_heat_demand_kWh',
-    'appliances_electrical_demand_kWh',
-    'lighting_electrical_demand_kWh'
-]
 
 # # Function to extract demand data from the GeoJSON
 # def extract_demand_data(geojson_data, demand_type, period_start, period_end):
@@ -122,6 +114,14 @@ demand_types = [
 #     plot_comparison(selected_buildings, baseline_summer, efficient_summer, range(period_start_summer, period_end_summer), f'One Week in Summer - Comparison - {demand_type}')
 #
 
+# Define the demand types to be plotted
+demand_types = [
+    'heating_demand_kWh',
+    'cooling_demand_kWh',
+    'domestic_hot_water_heat_demand_kWh',
+    'appliances_electrical_demand_kWh',
+    'lighting_electrical_demand_kWh'
+]
 
 # Function to extract and sum demand data from the GeoJSON
 def extract_and_sum_demand_data(geojson_data, demand_type, period_start, period_end):

hub/city_model_structure/building.py

@@ -810,37 +810,37 @@ 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],
+    # 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],
+    #                                          '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,
+    #                                          '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],
+    #                                                    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]},
+    #                                          'west': [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]},
-                                 cte.YEAR: {'north': [0],
+    #                              cte.YEAR: {'north': [0],
-                                            'east': [0],
+    #                                         'east': [0],
-                                            'south': [1.212544],
+    #                                         'south': [1.212544],
-                                            'west': [0]}
+    #                                         'west': [0]}
-                                 }
+    #                              }
-
+    #
-    # Add other systems whenever new ones appear
+    # # Add other systems whenever new ones appear
-    if self.energy_systems is None:
+    # if self.energy_systems is None:
-      return self._onsite_electrical_production
+    #   return self._onsite_electrical_production
-    for energy_system in self.energy_systems:
+    # for energy_system in self.energy_systems:
-      for generation_system in energy_system.generation_systems:
+    #   for generation_system in energy_system.generation_systems:
-        if generation_system.system_type == cte.PHOTOVOLTAIC:
+    #     if generation_system.system_type == cte.PHOTOVOLTAIC:
-          if generation_system.electricity_efficiency is not None:
+    #       if generation_system.electricity_efficiency is not None:
-            _efficiency = float(generation_system.electricity_efficiency)
+    #         _efficiency = float(generation_system.electricity_efficiency)
-          else:
+    #       else:
-            _efficiency = 0
+    #         _efficiency = 0
-          self._onsite_electrical_production = {}
+    #       self._onsite_electrical_production = {}
-          for _key in self.roofs[0].global_irradiance.keys():
+    #       for _key in self.roofs[0].global_irradiance.keys():
-            _results = [0 for _ in range(0, len(self.roofs[0].global_irradiance[_key]))]
+    #         _results = [0 for _ in range(0, len(self.roofs[0].global_irradiance[_key]))]
-            for surface in self.roofs:
+    #         for surface in self.roofs:
-              if _key in orientation_losses_factor:
+    #           if _key in orientation_losses_factor:
-                _results = [x + y * _efficiency * surface.perimeter_area
+    #             _results = [x + y * _efficiency * surface.perimeter_area
-                            * surface.solar_collectors_area_reduction_factor * z
+    #                         * surface.solar_collectors_area_reduction_factor * z
-                            for x, y, z in zip(_results, surface.global_irradiance[_key],
+    #                         for x, y, z in zip(_results, surface.global_irradiance[_key],
-                                               orientation_losses_factor[_key]['south'])]
+    #                                            orientation_losses_factor[_key]['south'])]
-            self._onsite_electrical_production[_key] = _results
+    #         self._onsite_electrical_production[_key] = _results
     return self._onsite_electrical_production
 
   @property

hub/data/energy_systems/montreal_future_systems.xml

@@ -1410,8 +1410,8 @@
 			<demand>cooling</demand>
 		</demands>
 		<components>
-			<generation_id>23</generation_id>
+			<generation_id>25</generation_id>
-			<generation_id>16</generation_id>
+			<generation_id>17</generation_id>
 		</components>
 	</system>
     <system>

main.py

@@ -14,13 +14,14 @@ import subprocess
 from hub.imports.results_factory import ResultFactory
 from hub.imports.energy_systems_factory import EnergySystemsFactory
 from scripts.energy_system_sizing_and_simulation_factory import EnergySystemsSimulationFactory
+from scripts.solar_angles import CitySolarAngles
 import hub.helpers.constants as cte
 from hub.exports.exports_factory import ExportsFactory
-from scripts.solar_angles import CitySolarAngles
+from scripts.pv_sizing_and_simulation import PVSizingSimulation
-from scripts.radiation_tilted import RadiationTilted
+import pandas as pd
-
 import geopandas as gpd
 import json
+
 #%% # -----------------------------------------------
 # Specify the GeoJSON file path
 #%% # -----------------------------------------------
@@ -41,7 +42,7 @@ cost_analysis_output_path.mkdir(parents=True, exist_ok=True)
 #%% # -----------------------------------------------
 
 geojson_file_path_baseline = output_path / 'updated_buildings_with_all_data_baseline.geojson'
-geojson_file_path_2024 = output_path / 'updated_buildings_with_all_data.geojson'
+geojson_file_path_2024 = output_path / 'updated_buildings_with_all_data_2024.geojson'
 with open(geojson_file_path_baseline , 'r') as f:
     building_type_data = json.load(f)
 with open(geojson_file_path_2024, 'r') as f:
@@ -54,47 +55,14 @@ city = GeometryFactory('geojson',
                        year_of_construction_field='year_built',
                        function_field='building_type',
                        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()
 # #energy plus is not going to be processed here, as demand has been obtained before
 # energy_plus_workflow(city)
-#%%-----------------------------------------------
-# """SRA"""
-#%% # -----------------------------------------------
-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()
-solar_angles = CitySolarAngles(city.name,
-                               city.latitude,
-                               city.longitude,
-                               tilt_angle=45,
-                               surface_azimuth_angle=180).calculate
-for building in city.buildings:
-  ghi = [x / cte.WATTS_HOUR_TO_JULES for x in building.roofs[0].global_irradiance[cte.HOUR]]
-  RadiationTilted(building,
-                  solar_angles,
-                  tilt_angle=45,
-                  ghi=ghi).enrich()
-# building_names = []
-# for building in city.buildings:
-#   building_names.append(building.name)
-#
-# df = pd.DataFrame(columns=building_names)
-# df1 = pd.DataFrame(columns=building_names)
-# print('test')
-# for building in city.buildings:
-#   # if building.name in selected_buildings_list:
-#   df[f'{building.name}'] = building.roofs[0].global_irradiance[cte.HOUR]
-#   df1[f'{building.name}'] = building.roofs[0].global_irradiance_tilted[cte.HOUR]
-#
-# df.to_csv('solar_radiation_horizontal_selected_buildings.csv')
-# df1.to_csv('solar_radiation_tilted_selected_buildings.csv')
-
 #%% # -----------------------------------------------
 #"""Enrich city with geojson file data"""
 #%% # -----------------------------------------------
@@ -169,50 +137,616 @@ print('test')
 #%% # -----------------------------------------------
 
 for building in city.buildings:
-  building.energy_systems_archetype_name = 'system 1 electricity pv'
+  building.energy_systems_archetype_name = 'system 1 electricity'
-
 
 EnergySystemsFactory('montreal_custom', city).enrich()
-# for building in city.buildings:
+
-#   energy_systems = building.energy_systems
+def baseline_to_dict(building):
-#   for energy_system in energy_systems:
-#     generation_units = energy_system.generation_systems
-#     if cte.HEATING in energy_system.demand_types:
-#       for generation_unit in generation_units:
-#         generation_unit.heat_efficiency = 0.96
-def to_dict(building, total_floor_area):
     return {
-        'roof_area': building.floor_area,
+        'heating_consumption_kWh': [x / (cte.WATTS_HOUR_TO_JULES * 1000) for x in building.heating_consumption[cte.HOUR]],
-        'total_floor_area': total_floor_area,
+        'cooling_consumption_kWh':[x / (cte.WATTS_HOUR_TO_JULES * 1000) for x in  building.cooling_consumption[cte.HOUR]],
-        'heating_consumption_kWh': [x / (cte.WATTS_HOUR_TO_JULES * 1000) for x in building.energy_consumption_breakdown[cte.HOUR]],
+        'domestic_hot_water_consumption_kWh': [x / (cte.WATTS_HOUR_TO_JULES * 1000) for x in building.domestic_hot_water_consumption[cte.HOUR]],
-        'cooling_consumption_kWh':[x / (cte.WATTS_HOUR_TO_JULES * 1000) for x in  building.cooling_demand[cte.HOUR]],
-        'domestic_hot_water_consumption_kWh': [x / (cte.WATTS_HOUR_TO_JULES * 1000) for x in building.domestic_hot_water_heat_demand[cte.HOUR]],
         'appliances_consumption_kWh':[x / (cte.WATTS_HOUR_TO_JULES * 1000) for x in  building.appliances_electrical_demand[cte.HOUR]],
         'lighting_consumption_kWh': [x / (cte.WATTS_HOUR_TO_JULES * 1000) for x in building.lighting_electrical_demand[cte.HOUR]]
     }
 buildings_dic={}
+for building in city.buildings:
+    buildings_dic[building.name]=baseline_to_dict(building)
+scenario={}
+
+scenario['baseline']=buildings_dic
+print("Scenario 1: Baseline is performed successfully")
+
+
+
+del city
+del buildings_dic
+del building_type_data
+
+#%%-----------------------------------------------
+# Scenario 2
+#%% # -----------------------------------------------
+
+# Create city object from GeoJSON file
+city = GeometryFactory('geojson',
+                       path=geojson_file_path_2024,
+                       height_field='maximum_roof_height',
+                       year_of_construction_field='year_built',
+                       function_field='building_type',
+                       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()
+
+enrich_buildings_with_geojson_data (building_type_data_2024, city)
+
+def to_dict(building,hourly_pv):
+    return {
+        'heating_consumption_kWh': [x / (cte.WATTS_HOUR_TO_JULES * 1000) for x in building.heating_consumption[cte.HOUR]],
+        'cooling_consumption_kWh':[x / (cte.WATTS_HOUR_TO_JULES * 1000) for x in  building.cooling_consumption[cte.HOUR]],
+        'domestic_hot_water_consumption_kWh': [x / (cte.WATTS_HOUR_TO_JULES * 1000) for x in building.domestic_hot_water_consumption[cte.HOUR]],
+        'appliances_consumption_kWh':[x / (cte.WATTS_HOUR_TO_JULES * 1000) for x in  building.appliances_electrical_demand[cte.HOUR]],
+        'lighting_consumption_kWh': [x / (cte.WATTS_HOUR_TO_JULES * 1000) for x in building.lighting_electrical_demand[cte.HOUR]],
+        'hourly_pv_kWh': [x /(cte.WATTS_HOUR_TO_JULES * 1000) for x in hourly_pv]
+    }
+buildings_dic={}
 
 
 for building in city.buildings:
+  building.energy_systems_archetype_name = 'system 1 electricity pv'
 
-    total_floor_area = 0
+EnergySystemsFactory('montreal_custom', city).enrich()
-    for thermal_zone in building.thermal_zones_from_internal_zones:
+# #%%-----------------------------------------------
-        total_floor_area += thermal_zone.total_floor_area
+# # """SRA"""
+# #%% # -----------------------------------------------
+ExportsFactory('sra', city, output_path).export()
+sra_path = (output_path / f'{city.name}_sra.xml').resolve()
+subprocess.run(['sra', str(sra_path)])
+ResultFactory('sra', city, output_path).enrich()
+solar_angles = CitySolarAngles(city.name,
+                               city.latitude,
+                               city.longitude,
+                               tilt_angle=45,
+                               surface_azimuth_angle=180).calculate
+df = pd.DataFrame()
+df.index = ['yearly lighting (kWh)', 'yearly appliance (kWh)', 'yearly heating (kWh)', 'yearly cooling (kWh)',
+            'yearly dhw (kWh)', 'roof area (m2)', 'used area for pv (m2)', 'number of panels', 'pv production (kWh)']
+for building in city.buildings:
+  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()
+  yearly_lighting = building.lighting_electrical_demand[cte.YEAR][0] / 1000
+  yearly_appliance = building.appliances_electrical_demand[cte.YEAR][0] / 1000
+  yearly_heating = building.heating_demand[cte.YEAR][0] / (3.6e6 * 3)
+  yearly_cooling = building.cooling_demand[cte.YEAR][0] / (3.6e6 * 4.5)
+  yearly_dhw = building.domestic_hot_water_heat_demand[cte.YEAR][0] / 1000
+  roof_area = building.roofs[0].perimeter_area
+  used_roof = pv_sizing_simulation.available_space()
+  number_of_pv_panels = pv_sizing_simulation.total_number_of_panels
+  yearly_pv = building.onsite_electrical_production[cte.YEAR][0] / (3.6e6)
+  hourly_pv = building.onsite_electrical_production[cte.HOUR]
+  df[f'{building.name}'] = [yearly_lighting, yearly_appliance, yearly_heating, yearly_cooling, yearly_dhw, roof_area,
+                            used_roof, number_of_pv_panels, yearly_pv]
+  buildings_dic[building.name]=to_dict(building,hourly_pv)
+
+# %%-----------------------------------------------
+# """South facing facades"""
+# %% # -----------------------------------------------
+# Function to convert radians to degrees
+import math
+def radians_to_degrees(radians):
+    return radians * (180 / math.pi)
+# Step 1: Create the walls_id dictionary
+walls_id={}
+
+for building in city.buildings:
+    ids = {}
+    for walls in building.walls:
+        id=walls.id
+        azimuth_degree=radians_to_degrees(float(walls.azimuth))
+        if azimuth_degree>90.0 or azimuth_degree <float(-90.0):
+            ids[id]= {
+                'azimuth': azimuth_degree,
+                'global_irradiance': walls.global_irradiance[cte.HOUR],
+                'area': walls.perimeter_area
+            }
+    walls_id[building.name] = ids
+
+# Step 2: Calculate pv_on_facade for each wall
+for building_id, ids in walls_id.items():
+    for wall_id, wall_data in ids.items():
+        if 'global_irradiance' in wall_data:
+            ghi = [x / cte.WATTS_HOUR_TO_JULES/1000 for x in wall_data['global_irradiance']]
+            wall_data['pv_on_facade'] = [x * 0.6 * wall_data['area']*0.22 for x in ghi]
 
 
 
 
-# for building in city.buildings:
+walls_dic = output_path / 'walls_id.json'
-#   print(building.heating_demand[cte.YEAR][0] / 3.6e6)
+with open(walls_dic , 'w') as json_file:
-#   print(building.name)
+    json.dump(walls_id, json_file, indent=4)
-#   total_floor_area = 0
+
-#   for thermal_zone in building.thermal_zones_from_internal_zones:
+import pandas as pd
-#     total_floor_area += thermal_zone.total_floor_area
+#### EXPORT
-#   print(building.heating_demand[cte.YEAR][0] / (3.6e6 * total_floor_area))
+# Convert walls_id dictionary to a DataFrame
+# Convert walls_id dictionary to DataFrames for static and hourly data
+# def convert_walls_id_to_dfs(walls_id):
+#     static_data = {}
+#     hourly_data = {}
+#
+#     for building_id, ids in walls_id.items():
+#         for wall_id, wall_data in ids.items():
+#             # Static data
+#             static_data[f"{building_id}_{wall_id}_azimuth"] = wall_data.get('azimuth', None)
+#             static_data[f"{building_id}_{wall_id}_area"] = wall_data.get('area', None)
+#
+#             if 'pv_on_facade' in wall_data:
+#                 hourly_data[f"{building_id}_{wall_id}_pv_on_facade"] = wall_data['pv_on_facade']
+#
+#     # Create DataFrames
+#     static_df = pd.DataFrame([static_data])
+#     hourly_df = pd.DataFrame(hourly_data)
+#
+#     return static_df, hourly_df
 
 
+# output_path_walls_id_dic =output_path / 'walls_id_data.xlsx'
+#
+# static_df, hourly_df = convert_walls_id_to_dfs(walls_id)
+# with pd.ExcelWriter(output_path_walls_id_dic) as writer:
+#     static_df.to_excel(writer, sheet_name='Static Data', index=False)
+#     hourly_df.to_excel(writer, sheet_name='Hourly Data', index=False)
+
+# print(f"Data successfully exported to {output_path}")
+# # Save the DataFrame to an Excel file
+
+
+
+df.to_csv(output_path / 'pv.csv')
+
+scenario['efficient with PV']=buildings_dic
+print("Scenario 2: efficient with PV run successfully")
+
+# #%%-----------------------------------------------
+# # Scenario 3
+# #%% # -----------------------------------------------
+#
 # for building in city.buildings:
 #   building.energy_systems_archetype_name = 'PV+4Pipe+DHW'
 # EnergySystemsFactory('montreal_future', city).enrich()
+# buildings_dic = {}
 # for building in city.buildings:
 #   EnergySystemsSimulationFactory('archetype13', building=building, output_path=simulation_results_path).enrich()
+#   buildings_dic[building.name] = to_dict(building, hourly_pv)
+# scenario['efficient with PV+4Pipe+DHW']=buildings_dic
+# print("Scenario 3: efficient with PV+4Pipe+DHW run successfully")
+#
+# def extract_HP_size(building):
+#     dic={
+#         # Heat Pump Rated Heating and Cooling Output
+#         'hp_heat_size': building.energy_systems[1].generation_systems[1].nominal_heat_output/1000,
+#         'hp_cooling_output':  building.energy_systems[1].generation_systems[1].nominal_cooling_output/1000,
+#         # Boiler Rated Heat Output
+#         'boiler_heat_output': building.energy_systems[1].generation_systems[0].nominal_heat_output/1000,
+#         # TES characteristics
+#         'tes_volume':building.energy_systems[1].generation_systems[0].energy_storage_systems[0].volume,
+#         'tes_height':building.energy_systems[1].generation_systems[0].energy_storage_systems[0].height,
+#         # DHW HP
+#         'dhw_hp_heat_output': building.energy_systems[-1].generation_systems[0].nominal_heat_output/1000,
+#         # DHW TES Characteristics
+#         'dhw_tes_volume':  building.energy_systems[-1].generation_systems[0].energy_storage_systems[0].volume,
+#         'dhw_tes_height':  building.energy_systems[-1].generation_systems[0].energy_storage_systems[0].height,
+#     }
+#
+#
+#     return dic
+# HPs={}
+# for building in city.buildings:
+#     HPs[building.name]=extract_HP_size(building)
+
+clusters=pd.read_csv(output_path/'clusters.csv')
+
+
+# Step 2: Extract the demand data for each building
+def extract_building_demand(city):
+    building_demand = {}
+    for building in city.buildings:
+        demands = {
+            'heating_demand': [x / (1000 * cte.WATTS_HOUR_TO_JULES) for x in building.heating_demand[cte.HOUR]],
+            'cooling_demand': [x / (1000 * cte.WATTS_HOUR_TO_JULES) for x in building.cooling_demand[cte.HOUR]],
+            'domestic_hot_water_demand': [x / (1000 * cte.WATTS_HOUR_TO_JULES) for x in building.domestic_hot_water_heat_demand[cte.HOUR]],
+            'appliances_electrical_demand': [x / (1000 * cte.WATTS_HOUR_TO_JULES) for x in building.appliances_electrical_demand[cte.HOUR]],
+            'lighting_electrical_demand': [x / (1000 * cte.WATTS_HOUR_TO_JULES) for x in building.lighting_electrical_demand[cte.HOUR]]
+        }
+        building_demand[building.name] = demands
+    return building_demand
+
+# Step 3: Sum the demand types for each cluster
+def sum_demands_by_cluster(building_demand, clusters, demand_types):
+    cluster_demands = {cluster: {demand_type: [0] * 8760 for demand_type in demand_types} for cluster in clusters['cluster'].unique()}
+
+    for _, row in clusters.iterrows():
+        building_id = str(row['id'])
+        cluster = row['cluster']
+        if building_id in building_demand:
+            for demand_type in demand_types:
+                cluster_demands[cluster][demand_type] = [sum(x) for x in zip(cluster_demands[cluster][demand_type], building_demand[building_id][demand_type])]
+
+    return cluster_demands
+
+
+def plot_demands_by_cluster(cluster_demands, demand_types, output_folder):
+    import os
+    os.makedirs(output_folder, exist_ok=True)
+
+    for cluster, demands in cluster_demands.items():
+        plt.figure(figsize=(15, 10))
+        for demand_type in demand_types:
+            plt.plot(demands[demand_type], label=demand_type)
+
+        plt.title(f'Summed Demands for Cluster {cluster}')
+        plt.xlabel('Hour of the Year')
+        plt.ylabel('Demand (kWh)')
+        plt.legend(loc='upper right')
+        plt.grid(True)
+        plt.tight_layout()
+        plt.savefig(os.path.join(output_folder, f'cluster_{cluster}_summed_demands.png'))
+        plt.close()
+
+
+# Example usage
+demand_types = [
+    'heating_demand',
+    'cooling_demand',
+    'domestic_hot_water_demand',
+    'appliances_electrical_demand',
+    'lighting_electrical_demand'
+]
+
+# Extract the building demand data
+building_demand = extract_building_demand(city)
+cluster_demands = sum_demands_by_cluster(building_demand, clusters, demand_types)
+# Create a DataFrame to export the results
+cluster_demands_df = {f"{cluster}_{demand_type}": data for cluster, demands in cluster_demands.items() for
+                      demand_type, data in demands.items()}
+cluster_demands_df = pd.DataFrame(cluster_demands_df)
+
+# Save the results to an Excel file
+
+cluster_demands_df.to_excel(output_path/'cluster_demands.xlsx', index=False)
+
+print(f"Clustered demand data successfully exported to {output_path}")
+
+
+
+#%%-----------------------------------------------
+# Scenario 4
+#%% # -----------------------------------------------
+
+del city
+del buildings_dic
+
+
+geojson_file_path_clusters= output_path / 'new.geojson'
+
+with open(geojson_file_path_clusters , 'r') as f:
+    building_type_data_new = json.load(f)
+
+# Create city object from GeoJSON file
+city = GeometryFactory('geojson',
+                       path=geojson_file_path_clusters,
+                       height_field='maximum_roof_height',
+                       year_of_construction_field='year_built',
+                       function_field='building_type',
+                       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()
+
+buildings_clusters={
+    1651: 4,
+    1662: 0,
+    1667: 1,
+    1674: 2,
+    1688: 3
+}
+
+for building_id in buildings_clusters:
+    cluster=buildings_clusters[building_id]
+    for idx, feature in enumerate(building_type_data_new['features']):
+        if feature['properties']['id'] == str(building_id):
+            building_type_data_new['features'][idx]['properties']['heating_demand_kWh']=cluster_demands[cluster]['heating_demand']
+            building_type_data_new['features'][idx]['properties']['cooling_demand_kWh'] = cluster_demands[cluster]['cooling_demand']
+            building_type_data_new['features'][idx]['properties']['domestic_hot_water_heat_demand_kWh'] = cluster_demands[cluster]['domestic_hot_water_demand']
+            building_type_data_new['features'][idx]['properties']['appliances_electrical_demand_kWh'] = cluster_demands[cluster]['appliances_electrical_demand']
+            building_type_data_new['features'][idx]['properties']['lighting_electrical_demand_kWh'] = cluster_demands[cluster]['lighting_electrical_demand']
+
+enrich_buildings_with_geojson_data (building_type_data_new, city)
+
+
+for building in city.buildings:
+  building.energy_systems_archetype_name = 'PV+4Pipe+DHW'
+EnergySystemsFactory('montreal_future', city).enrich()
+buildings_dic = {}
+for building in city.buildings:
+  EnergySystemsSimulationFactory('archetype13', building=building, output_path=simulation_results_path).enrich()
+  buildings_dic[building.name] = to_dict(building, hourly_pv)
+scenario['efficient with PV+4Pipe+DHW']=buildings_dic
+print("Scenario 4: efficient with PV+4Pipe+DHW run successfully for Clusters")
+
+def extract_HP_size(building):
+    dic={
+        # Heat Pump Rated Heating and Cooling Output
+        'hp_heat_size': building.energy_systems[1].generation_systems[1].nominal_heat_output/1000,
+        'hp_cooling_output':  building.energy_systems[1].generation_systems[1].nominal_cooling_output/1000,
+        # Boiler Rated Heat Output
+        'boiler_heat_output': building.energy_systems[1].generation_systems[0].nominal_heat_output/1000,
+        # TES characteristics
+        'tes_volume':building.energy_systems[1].generation_systems[0].energy_storage_systems[0].volume,
+        'tes_height':building.energy_systems[1].generation_systems[0].energy_storage_systems[0].height,
+        # DHW HP
+        'dhw_hp_heat_output': building.energy_systems[-1].generation_systems[0].nominal_heat_output/1000,
+        # DHW TES Characteristics
+        'dhw_tes_volume':  building.energy_systems[-1].generation_systems[0].energy_storage_systems[0].volume,
+        'dhw_tes_height':  building.energy_systems[-1].generation_systems[0].energy_storage_systems[0].height,
+    }
+
+
+    return dic
+HPs={}
+for building in city.buildings:
+    HPs[building.name]=extract_HP_size(building)
+
+#%%-------------------------------------------------------
+#""""EXPORTERS"""
+#%%-------------------------------------------------------
+
+
+# Convert the dictionary to a DataFrame
+df = pd.DataFrame.from_dict(HPs, orient='index')
+
+# Save the DataFrame to an Excel file
+output_path_HPs =output_path/ 'HPs_data_sc4.xlsx'
+df.to_excel(output_path_HPs, index_label='building_id')
+
+print(f"Data successfully exported to {output_path}")
+
+
+
+#%%-------------------------------------------------------
+#""""EXPORTERS"""
+#%%-------------------------------------------------------
+
+
+# Convert the dictionary to a DataFrame
+df = pd.DataFrame.from_dict(HPs, orient='index')
+
+# Save the DataFrame to an Excel file
+output_path_HPs =output_path/ 'HPs_data.xlsx'
+df.to_excel(output_path_HPs, index_label='building_id')
+
+print(f"Data successfully exported to {output_path}")
+
+import pandas as pd
+
+districts_demands={}
+def extract_and_sum_demand_data(scenario, demand_types):
+
+
+    # Conversion factor constant
+    conversion_factor = 1 / (cte.WATTS_HOUR_TO_JULES * 1000)
+
+    # Loop through each scenario
+    for scenario_key, buildings in scenario.items():
+        # Loop through each building in the scenario
+        # Initialize an empty dictionary to store the district demand sums
+        district_demand = {demand_type: [0] * 8760 for demand_type in demand_types}
+        district_demand['hourly_pv_kWh']= [0] * 8760
+        for building_id, building_data in buildings.items():
+            # Loop through each demand type and sum up the data
+            for demand_type in demand_types:
+                if demand_type in building_data:
+                    district_demand[demand_type] = [sum(x) for x in zip(district_demand[demand_type], building_data[demand_type])]
+
+                # If PV data is available and relevant
+                if scenario_key == "efficient with PV":
+                    district_demand['hourly_pv_kWh'] = [sum(x) for x in zip(district_demand['hourly_pv_kWh'], building_data['hourly_pv_kWh'])]
+            if scenario_key == 'efficient with PV+4Pipe+DHW':
+                district_demand['hourly_pv_kWh'] = districts_demands["efficient with PV"]['hourly_pv_kWh']
+            districts_demands[scenario_key]=district_demand
+
+    return districts_demands
+
+# Example usage
+# Assuming 'scenario' is a dictionary with the required structure and 'cte' is defined somewhere with WATTS_HOUR_TO_JULES constant
+demand_types = [
+    'heating_consumption_kWh',
+    'cooling_consumption_kWh',
+    'domestic_hot_water_consumption_kWh',
+    'appliances_consumption_kWh',
+    'lighting_consumption_kWh',
+    # 'hourly_pv_kWh'  # Include this only if you want to consider PV data
+]
+
+# # Call the function with your scenario data
+district_demand = extract_and_sum_demand_data(scenario, demand_types)
+#
+# """"EXPORTERS"""
+# import pandas as pd
+#
+#
+# Export the DataFrame to an Excel file
+excel_file_path = r'C:\Users\a_gabald\PycharmProjects\summer_course_2024\out_files\districts_balance.xlsx'
+# df.to_excel(excel_file_path, index=True, index_label='Building')
+
+# Create an Excel writer object
+with pd.ExcelWriter(excel_file_path, engine='xlsxwriter') as writer:
+    for scenarios,demands in district_demand.items():
+        # Convert demands to a DataFrame
+        df_demands = pd.DataFrame(demands)
+        # Convert building_id to string and check its length
+        sheet_name = str(scenarios)
+        if len(sheet_name) > 31:
+            sheet_name = sheet_name[:31]  # Truncate to 31 characters if necessary
+        # Write the DataFrame to a specific sheet named after the building_id
+        df_demands.to_excel(writer, sheet_name=sheet_name, index=False)
+
+
+print("district balance data is exported successfully")
+
+
+import pandas as pd
+
+# Assuming your scenario dictionary is already defined as follows:
+# scenario = {
+#     'baseline': { ... },
+#     'efficient with PV': { ... }
+# }
+
+
+def dict_to_df_col_wise(building_data):
+    """
+    Converts a dictionary of building data to a DataFrame.
+
+    Args:
+        building_data (dict): Dictionary containing building data where keys are building ids and values are dictionaries
+                              with hourly data for various demand types.
+
+    Returns:
+        pd.DataFrame: DataFrame with columns for each building and demand type.
+    """
+    # Create a dictionary to hold DataFrames for each demand type
+    df_dict= {}
+
+    # Loop over each building
+    for building_id, data in building_data.items():
+        # Create a DataFrame for this building's data
+        building_df = pd.DataFrame(data)
+
+        # Rename columns to include building_id
+        building_df.columns = [f"{building_id}_{col}" for col in building_df.columns]
+
+        # Add this DataFrame to the dictionary
+        df_dict[building_id] = building_df
+
+    # Concatenate all building DataFrames column-wise
+    result_df = pd.concat(df_dict.values(), axis=1)
+
+    return result_df
+
+# Create DataFrames for each scenario
+baseline_df = dict_to_df_col_wise(scenario['baseline'])
+efficient_with_pv_df = dict_to_df_col_wise(scenario['efficient with PV'])
+efficient_with_pv_hps = dict_to_df_col_wise(scenario['efficient with PV+4Pipe+DHW'])
+
+
+# Write the DataFrames to an Excel file with two separate sheets
+with pd.ExcelWriter(r'C:\Users\a_gabald\PycharmProjects\summer_course_2024\out_files\scenario_data.xlsx') as writer:
+    baseline_df.to_excel(writer, sheet_name='baseline', index=True)
+    efficient_with_pv_df.to_excel(writer, sheet_name='efficient with PV', index=True)
+    efficient_with_pv_hps.to_excel(writer, sheet_name='efficient with HPs_2', index=True)
+
+print("hourly data has been successfully exported per building to scenario_data.xlsx")
+
+
+
+import pandas as pd
+import matplotlib.pyplot as plt
+import numpy as np
+
+
+def convert_hourly_to_monthly(hourly_data):
+    """
+    Converts hourly data to monthly data by summing up the values for each month.
+
+    Args:
+        hourly_data (list): List of hourly data (length 8760).
+
+    Returns:
+        list: List of monthly data (length 12).
+    """
+    hourly_series = pd.Series(hourly_data, index=pd.date_range(start='1/1/2023', periods=8760, freq='H'))
+    monthly_data = hourly_series.resample('M').sum()
+    return monthly_data.tolist()
+
+import os
+def plot_stacked_demands_vs_pv(district_demand, demand_types, output_path, pv_type='hourly_pv_kWh'):
+    """
+    Plots the stacked monthly demand for each scenario and compares it to the PV data.
+
+    Args:
+        district_demand (dict): Dictionary with scenario keys and demand data.
+        demand_types (list): List of demand types to plot.
+        output_path (str): Path to save the plots.
+        pv_type (str): The PV data type to compare against.
+    """
+    os.makedirs(output_path, exist_ok=True)
+
+    for scenario_key, demand_data in district_demand.items():
+        # Convert hourly data to monthly data for each demand type
+        monthly_data = {demand_type: convert_hourly_to_monthly(demand_data[demand_type]) for demand_type in
+                        demand_types}
+        monthly_pv = convert_hourly_to_monthly(demand_data.get(pv_type, [0] * 8760))
+
+        # Create a DataFrame for easier plotting
+        combined_data = pd.DataFrame(monthly_data)
+        combined_data['Month'] = range(1, 13)
+        combined_data['PV'] = monthly_pv
+
+        # Plotting
+        fig, ax1 = plt.subplots(figsize=(14, 8))
+
+        # Plot stacked demands
+        combined_data.set_index('Month', inplace=True)
+        combined_data[demand_types].plot(kind='bar', stacked=True, ax=ax1, colormap='tab20')
+
+        ax1.set_xlabel('Month')
+        ax1.set_ylabel('Energy Demand (kWh)')
+        ax1.set_title(f'Monthly Energy Demand and PV Generation for {scenario_key}')
+
+        # Plot PV data on the secondary y-axis
+        ax2 = ax1.twinx()
+        ax2.plot(combined_data.index, combined_data['PV'], color='black', linestyle='-', marker='o',
+                 label='PV Generation')
+        ax2.set_ylabel('PV Generation (kWh)')
+
+        # Add legends
+        ax1.legend(loc='upper left')
+        ax2.legend(loc='upper right')
+
+        ax1.set_xticks(combined_data.index)
+        ax1.set_xticklabels(['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'])
+
+        # Save the plot
+        plt.savefig(os.path.join(output_path, f'{scenario_key}_monthly_demand_vs_pv.png'))
+        plt.close()
+
+
+# Example usage
+# district_demand = extract_and_sum_demand_data(scenario, demand_types)
+
+# Specify the demand types and PV type
+demand_types = [
+    'heating_consumption_kWh',
+    'cooling_consumption_kWh',
+    'domestic_hot_water_consumption_kWh',
+    'appliances_consumption_kWh',
+    'lighting_consumption_kWh'
+]
+
+# Plot the data
+plot_stacked_demands_vs_pv(district_demand, demand_types, output_path)
+# Plot the data
+print('test')
+import csv

processor.py

@@ -3,16 +3,18 @@ from pathlib import Path
 import pandas as pd
 import geopandas as gpd
 
-
+#%%-----------------------------------------------
+# """This code takes different results from energy plus with different building uses and obtains a unique geojson that combines all of them"""
+#%% # -----------------------------------------------
 #define output path
 output_path = (Path(__file__).parent / 'out_files').resolve()
 output_path.mkdir(parents=True, exist_ok=True)
 # Define output folders
 output_folders = ['building_type', 'building_type_2', 'building_type_3']
 #output paths containing energy+ results are read
-path1=r'C:\Users\a_gabald\PycharmProjects\summer_course_2024\out_files\building_type\updated_buildings.geojson'
+path1=output_path / '1990_and_2019'/'building_type' / 'updated_buildings.geojson'
-path2=r'C:\Users\a_gabald\PycharmProjects\summer_course_2024\out_files\building_type_2\updated_buildings.geojson'
+path2=output_path  / '1990_and_2019'/ 'building_type_2' / 'updated_buildings.geojson'
-path3=r'C:\Users\a_gabald\PycharmProjects\summer_course_2024\out_files\building_type_3\updated_buildings.geojson'
+path3=output_path  / '1990_and_2019'/'building_type_3' / 'updated_buildings.geojson'
 with open(path1, 'r') as f:
     building_type_data = json.load(f)
 with open(path2, 'r') as f:

scripts/energy_system_sizing_and_simulation_factory.py

@@ -4,8 +4,8 @@ SPDX - License - Identifier: LGPL - 3.0 - or -later
 Copyright © 2022 Concordia CERC group
 Project Coder Saeed Ranjbar saeed.ranjbar@mail.concordia.ca
 """
-
+from scripts.system_simulation_models.archetype13_electric import Archetype13
-from scripts.system_simulation_models.archetype13 import Archetype13
+# 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

scripts/system_simulation_models/archetype1.py

@@ -16,7 +16,8 @@ class Archetype1:
     self._domestic_hot_water_peak_load = building.domestic_hot_water_peak_load[cte.YEAR][0]
     self._hourly_heating_demand = [0] + [demand / 3600 for demand in building.heating_demand[cte.HOUR]]
     self._hourly_cooling_demand = [demand / 3600 for demand in building.cooling_demand[cte.HOUR]]
-    self._hourly_dhw_demand = building.domestic_hot_water_heat_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 = [0] + building.external_temperature[cte.HOUR]
     self.results = {}
@@ -27,9 +28,9 @@ class Archetype1:
     heat_pump = self._hvac_system.generation_systems[0]
     boiler = self._hvac_system.generation_systems[1]
     thermal_storage = heat_pump.energy_storage_systems[0]
-    heat_pump.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 / 3600)
+    heat_pump.nominal_cooling_output = round(self._cooling_peak_load)
-    boiler.nominal_heat_output = round(0.5 * self._heating_peak_load / 3600)
+    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))
@@ -60,7 +61,7 @@ class Archetype1:
     (t_sup_hp, t_tank, t_ret, m_ch, m_dis, q_hp, q_boiler, hp_cop,
      hp_electricity, boiler_gas, boiler_consumption, t_sup_boiler, heating_consumption) = [variables[name] for name in
                                                                                            variable_names]
-    t_tank[0] = 55
+    t_tank[0] = 30
     dt = 3600
     hp_heating_cap = hp.nominal_heat_output
     hp_efficiency = float(hp.heat_efficiency)

scripts/system_simulation_models/archetype13_electric.py

@@ -22,12 +22,13 @@ class Archetype13:
     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._t_out = [15]*8760
+    # building.external_temperature[cte.HOUR]
     self.results = {}
     self.dt = 900
 
-  def hvac_sizing(self):
+  def hvac_sizing(self,storage_factor=3):
-    storage_factor = 3
+    storage_factor = storage_factor
     heat_pump = self._hvac_system.generation_systems[1]
     boiler = self._hvac_system.generation_systems[0]
     thermal_storage = boiler.energy_storage_systems[0]
@@ -39,8 +40,8 @@ class Archetype13:
       (cte.WATER_HEAT_CAPACITY * cte.WATER_DENSITY * 25))
     return heat_pump, boiler, thermal_storage
 
-  def dhw_sizing(self):
+  def dhw_sizing(self,storage_factor=3):
-    storage_factor = 3
+    storage_factor = storage_factor
     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
@@ -51,7 +52,7 @@ class Archetype13:
 
   def heating_system_simulation(self):
     hp, boiler, tes = self.hvac_sizing()
-    cop_curve_coefficients = [float(coefficient) for coefficient in hp.heat_efficiency_curve.coefficients]
+    # 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)]
@@ -102,12 +103,13 @@ class Archetype13:
       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] +
+        hp_cop[i + 1] =4.07
-                         cop_curve_coefficients[1] * t_tank_fahrenheit +
+          # (1 / (cop_curve_coefficients[0] +
-                         cop_curve_coefficients[2] * t_tank_fahrenheit ** 2 +
+          #                cop_curve_coefficients[1] * t_tank_fahrenheit +
-                         cop_curve_coefficients[3] * t_out_fahrenheit +
+          #                cop_curve_coefficients[2] * t_tank_fahrenheit ** 2 +
-                         cop_curve_coefficients[4] * t_out_fahrenheit ** 2 +
+          #                cop_curve_coefficients[3] * t_out_fahrenheit +
-                         cop_curve_coefficients[5] * t_tank_fahrenheit * t_out_fahrenheit)) * hp_efficiency
+          #                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
@@ -119,7 +121,7 @@ class Archetype13:
         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)
+        # 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:
@@ -178,7 +180,7 @@ class Archetype13:
 
   def cooling_system_simulation(self):
     hp = self.hvac_sizing()[0]
-    eer_curve_coefficients = [float(coefficient) for coefficient in hp.cooling_efficiency_curve.coefficients]
+    # 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)]
@@ -216,12 +218,13 @@ class Archetype13:
       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] +
+        hp_cop[i] = 4
-                     eer_curve_coefficients[1] * t_sup_hp_fahrenheit +
+        # (1 / (eer_curve_coefficients[0] +
-                     eer_curve_coefficients[2] * t_sup_hp_fahrenheit ** 2 +
+        #              eer_curve_coefficients[1] * t_sup_hp_fahrenheit +
-                     eer_curve_coefficients[3] * t_out_fahrenheit +
+        #              eer_curve_coefficients[2] * t_sup_hp_fahrenheit ** 2 +
-                     eer_curve_coefficients[4] * t_out_fahrenheit ** 2 +
+        #              eer_curve_coefficients[3] * t_out_fahrenheit +
-                     eer_curve_coefficients[5] * t_sup_hp_fahrenheit * t_out_fahrenheit)) * cooling_efficiency / 3.41
+        #              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
@@ -251,7 +254,7 @@ class Archetype13:
 
   def dhw_system_simulation(self):
     hp, tes = self.dhw_sizing()
-    cop_curve_coefficients = [float(coefficient) for coefficient in hp.heat_efficiency_curve.coefficients]
+    # 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)]
@@ -298,12 +301,13 @@ class Archetype13:
       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] +
+        hp_cop[i] = 3
-                     cop_curve_coefficients[1] * t_out[i] +
+                     # (cop_curve_coefficients[0] +
-                     cop_curve_coefficients[2] * t_out[i] ** 2 +
+                     # cop_curve_coefficients[1] * t_out[i] +
-                     cop_curve_coefficients[3] * t_tank[i] +
+                     # cop_curve_coefficients[2] * t_out[i] ** 2 +
-                     cop_curve_coefficients[4] * t_tank[i] ** 2 +
+                     # cop_curve_coefficients[3] * t_tank[i] +
-                     cop_curve_coefficients[5] * t_tank[i] * t_out[i]) * float(hp.heat_efficiency)
+                     # 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