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/data/energy_systems/montreal_future_systems.xml

@@ -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>
 

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/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

@@ -29,7 +29,9 @@ residential_systems_percentage = {'system 1 gas': 100,
                                   'system 8 electricity': 0}
 
 residential_new_systems_percentage = {'PV+ASHP+GasBoiler+TES': 0,
-                                      'PV+4Pipe+DHW': 100,
+                                      'PV+4Pipe+DHW': 0,
+                                      'Central Heating+Unitary Cooling+Unitary DHW': 50,
+                                      'Central Heating+Unitary Cooling+Unitary DHW+PV': 50,
                                       'PV+ASHP+ElectricBoiler+TES': 0,
                                       'PV+GSHP+GasBoiler+TES': 0,
                                       'PV+GSHP+ElectricBoiler+TES': 0,

scripts/system_simulation_models/archetype13.py

@@ -19,7 +19,8 @@ class Archetype13:
     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_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):

scripts/system_simulation_models/archetypes14_15.py

@@ -0,0 +1,402 @@
+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()
+    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_sup_hp_fahrenheit = 1.8 * t_sup_hp[i + 1] + 32
+      t_out_fahrenheit = 1.8 * t_out[i + 1] + 32
+      if q_hp[i + 1] > 0:
+        hp_cop[i + 1] = (cop_curve_coefficients[0] +
+                         cop_curve_coefficients[1] * t_sup_hp_fahrenheit +
+                         cop_curve_coefficients[2] * t_sup_hp_fahrenheit ** 2 +
+                         cop_curve_coefficients[3] * t_out_fahrenheit +
+                         cop_curve_coefficients[4] * t_out_fahrenheit ** 2 +
+                         cop_curve_coefficients[5] * t_sup_hp_fahrenheit * t_out_fahrenheit)
+        hp_electricity[i + 1] = q_hp[i + 1] / 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()
+    eer_curve_coefficients = [float(coefficient) for coefficient in hp.cooling_efficiency_curve.coefficients]
+    cooling_efficiency = float(hp.cooling_efficiency)
+    number_of_ts = int(cte.HOUR_TO_SECONDS / self.dt)
+    demand = [0] + [x for x in self._hourly_cooling_demand for _ in range(number_of_ts)]
+    t_out = [0] + [x for x in self._t_out for _ in range(number_of_ts)]
+    hp.source_temperature = self._t_out
+    variable_names = ["t_sup_hp", "t_ret", "m", "q_hp", "hp_electricity", "hp_eer"]
+    num_hours = len(demand)
+    variables = {name: [0] * num_hours for name in variable_names}
+    (t_sup_hp, t_ret, m, q_hp, hp_electricity, hp_eer) = [variables[name] for name in variable_names]
+    t_ret[0] = 13
+
+    for i in range(1, len(demand)):
+      if demand[i] > 0:
+        m[i] = self._cooling_peak_load / (cte.WATER_HEAT_CAPACITY * 5 * cte.HOUR_TO_SECONDS)
+        if t_ret[i - 1] >= 13:
+          if demand[i] < 0.25 * self._cooling_peak_load / cte.HOUR_TO_SECONDS:
+            q_hp[i] = 0.25 * hp.nominal_cooling_output
+          elif demand[i] < 0.5 * self._cooling_peak_load / cte.HOUR_TO_SECONDS:
+            q_hp[i] = 0.5 * hp.nominal_cooling_output
+          else:
+            q_hp[i] = hp.nominal_cooling_output
+          t_sup_hp[i] = t_ret[i - 1] - q_hp[i] / (m[i] * cte.WATER_HEAT_CAPACITY)
+        else:
+          q_hp[i] = 0
+          t_sup_hp[i] = t_ret[i - 1]
+        if m[i] == 0:
+          t_ret[i] = t_sup_hp[i]
+        else:
+          t_ret[i] = t_sup_hp[i] + demand[i] / (m[i] * cte.WATER_HEAT_CAPACITY)
+      else:
+        m[i] = 0
+        q_hp[i] = 0
+        t_sup_hp[i] = t_ret[i -1]
+        t_ret[i] = t_ret[i - 1]
+      t_sup_hp_fahrenheit = 1.8 * t_sup_hp[i] + 32
+      t_out_fahrenheit = 1.8 * t_out[i] + 32
+      if q_hp[i] > 0:
+        hp_eer[i] = (eer_curve_coefficients[0] +
+                     eer_curve_coefficients[1] * t_sup_hp_fahrenheit +
+                     eer_curve_coefficients[2] * t_sup_hp_fahrenheit ** 2 +
+                     eer_curve_coefficients[3] * t_out_fahrenheit +
+                     eer_curve_coefficients[4] * t_out_fahrenheit ** 2 +
+                     eer_curve_coefficients[5] * t_sup_hp_fahrenheit * t_out_fahrenheit)
+        hp_electricity[i] = q_hp[i] / hp_eer[i]
+      else:
+        hp_eer[i] = 0
+        hp_electricity[i] = 0
+    hp_electricity_j = [(x * cte.WATTS_HOUR_TO_JULES) / number_of_ts for x in hp_electricity]
+    hp_hourly = []
+    hp_sum = 0
+    for i in range(1, len(demand)):
+      hp_sum += hp_electricity_j[i]
+      if (i - 1) % number_of_ts == 0:
+        hp_hourly.append(hp_sum)
+        hp_sum = 0
+    hp.energy_consumption[cte.COOLING] = {}
+    hp.energy_consumption[cte.COOLING][cte.HOUR] = hp_hourly
+    hp.energy_consumption[cte.COOLING][cte.MONTH] = MonthlyValues.get_total_month(
+      hp.energy_consumption[cte.COOLING][cte.HOUR])
+    hp.energy_consumption[cte.COOLING][cte.YEAR] = [
+      sum(hp.energy_consumption[cte.COOLING][cte.MONTH])]
+    self.results['Cooling Demand (W)'] = demand
+    self.results['HP Cooling Output (W)'] = q_hp
+    self.results['HP Cooling Supply Temperature'] = t_sup_hp
+    self.results['HP Cooling COP'] = hp_eer
+    self.results['HP Electricity Consumption'] = hp_electricity
+    self.results['Cooling Loop Flow Rate (kg/s)'] = m
+    self.results['Cooling Loop Return Temperature'] = t_ret
+    return hp_hourly
+
+  def dhw_system_simulation(self):
+    hp, tes = self.dhw_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] < 65:
+        q_hp[i] = hp_heating_cap
+      delta_t_hp = q_hp[i] * (self.dt / (cte.WATER_DENSITY * cte.WATER_HEAT_CAPACITY * v))
+      if demand[i] > 0:
+        dhw_needed = (demand[i] * cte.HOUR_TO_SECONDS) / (cte.WATER_HEAT_CAPACITY * t_tank[i] * cte.WATER_DENSITY)
+        m_dis[i] = dhw_needed * cte.WATER_DENSITY / cte.HOUR_TO_SECONDS
+        m_refill[i] = m_dis[i]
+      delta_t_freshwater = m_refill[i] * (t_tank[i] - freshwater_temperature) * (self.dt / (v * cte.WATER_DENSITY))
+      diff = delta_t_freshwater + delta_t_demand - delta_t_hp
+      if diff > 0:
+          if diff > 0:
+            power = diff * (cte.WATER_DENSITY * cte.WATER_HEAT_CAPACITY * v) / self.dt
+            if power <= float(tes.heating_coil_capacity):
+              q_coil[i] = power
+            else:
+              q_coil[i] = float(tes.heating_coil_capacity)
+      delta_t_coil = q_coil[i] * (self.dt / (cte.WATER_DENSITY * cte.WATER_HEAT_CAPACITY * v))
+
+      if q_hp[i] > 0:
+        m_ch[i] = q_hp[i] / (cte.WATER_HEAT_CAPACITY * hp_delta_t)
+        t_sup_hp[i] = (q_hp[i] / (m_ch[i] * cte.WATER_HEAT_CAPACITY)) + t_tank[i]
+      else:
+        m_ch[i] = 0
+        t_sup_hp[i] = t_tank[i]
+      t_sup_hp_fahrenheit = 1.8 * t_sup_hp[i] + 32
+      t_out_fahrenheit = 1.8 * t_out[i] + 32
+      if q_hp[i] > 0:
+        hp_cop[i] = (cop_curve_coefficients[0] +
+                     cop_curve_coefficients[1] * t_sup_hp_fahrenheit +
+                     cop_curve_coefficients[2] * t_sup_hp_fahrenheit ** 2 +
+                     cop_curve_coefficients[3] * t_out_fahrenheit +
+                     cop_curve_coefficients[4] * t_out_fahrenheit ** 2 +
+                     cop_curve_coefficients[5] * t_sup_hp_fahrenheit * t_out_fahrenheit)
+        hp_electricity[i] = q_hp[i] / 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()))

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):