energy_system_modelling_workflow/scripts/system_simulation.py

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)
        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) / 50e6
      boiler_consumption[i + 1] = q_aux[i + 1] / boiler_efficiency
      heating_consumption[i + 1] = boiler_consumption[i + 1] + hp_electricity[i + 1]
    data = list(zip(T, T_sup, T_ret, m_ch, m_dis, q_hp, hp_electricity, aux_fuel, q_aux, self.heating_demand))
    file_name = f'simulation_results_{self.name}.csv'
    with open(out_path / file_name, 'w', newline='') as csvfile:
      output_file = csv.writer(csvfile)
      # Write header
      output_file.writerow(['T', 'T_sup', 'T_ret', 'm_ch', 'm_dis', 'q_hp', 'hp_electricity', 'aux_fuel', 'q_aux', 'heating_demand'])
      # Write data
      output_file.writerows(data)
    return heating_consumption, hp_electricity, boiler_consumption, 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_consumption_disaggregated = disaggregated_consumption
    return self.building