s_ranjbar/energy_system_modelling_workflow
2
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge pull request 'fix: bugs in costing solved' (#6) from radiation_tilted into main

Browse Source 
Reviewed-on: https://nextgenerations-cities.encs.concordia.ca/gitea/s_ranjbar/energy_system_modelling_workflow/pulls/6
This commit is contained in:
Saeed Ranjbar 2024-05-28 10:37:51 -04:00
commit c14c88005d
14 changed files with 370 additions and 19 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

17
hub/city_model_structure/building_demand/surface.py
View File

@ -46,6 +46,7 @@ class Surface:
self._vegetation = None
self._percentage_shared = None
self._solar_collectors_area_reduction_factor = None
self._global_irradiance_tilted = {}
@property
def name(self):
@ -384,3 +385,19 @@ class Surface:
:param value: float
"""
self._solar_collectors_area_reduction_factor = value
@property
def global_irradiance_tilted(self) -> dict:
"""
Get global irradiance on a tilted surface in J/m2
:return: dict
"""
return self._global_irradiance_tilted
@global_irradiance_tilted.setter
def global_irradiance_tilted(self, value):
"""
Set global irradiance on a tilted surface in J/m2
:param value: dict
"""
self._global_irradiance_tilted = value

29
hub/city_model_structure/city_object.py
View File

@ -41,9 +41,10 @@ class CityObject:
self._ground_temperature = {}
self._global_horizontal = {}
self._diffuse = {}
self._beam = {}
self._direct_normal = {}
self._sensors = []
self._neighbours = None
self._beam = {}
@property
def level_of_detail(self) -> LevelOfDetail:
@ -238,20 +239,20 @@ class CityObject:
self._diffuse = value
@property
def beam(self) -> dict:
def direct_normal(self) -> dict:
"""
Get beam radiation surrounding the city object in J/m2
:return: dict{dict{[float]}}
"""
return self._beam
return self._direct_normal
@beam.setter
def beam(self, value):
@direct_normal.setter
def direct_normal(self, value):
"""
Set beam radiation surrounding the city object in J/m2
:param value: dict{dict{[float]}}
"""
self._beam = value
self._direct_normal = value
@property
def lower_corner(self):
@ -302,3 +303,19 @@ class CityObject:
Set the list of neighbour_objects and their properties associated to the current city_object
"""
self._neighbours = value
@property
def beam(self) -> dict:
"""
Get beam radiation surrounding the city object in J/m2
:return: dict{dict{[float]}}
"""
return self._beam
@beam.setter
def beam(self, value):
"""
Set beam radiation surrounding the city object in J/m2
:param value: dict{dict{[float]}}
"""
self._beam = value

36
hub/city_model_structure/energy_systems/pv_generation_system.py
View File

@ -26,6 +26,10 @@ class PvGenerationSystem(GenerationSystem):
self._width = None
self._height = None
self._electricity_power = None
self._tilt_angle = None
self._surface_azimuth = None
self._solar_altitude_angle = None
self._solar_azimuth_angle = None
@property
def nominal_electricity_output(self):
@ -202,3 +206,35 @@ class PvGenerationSystem(GenerationSystem):
:param value: float
"""
self._electricity_power = value
@property
def tilt_angle(self):
"""
Get tilt angle of PV system in degrees
:return: float
"""
return self._tilt_angle
@tilt_angle.setter
def tilt_angle(self, value):
"""
Set PV system tilt angle in degrees
:param value: float
"""
self._tilt_angle = value
@property
def surface_azimuth(self):
"""
Get surface azimuth angle of PV system in degrees. 0 is North
:return: float
"""
return self._surface_azimuth
@surface_azimuth.setter
def surface_azimuth(self, value):
"""
Set PV system tilt angle in degrees
:param value: float
"""
self._surface_azimuth = value

2
hub/exports/formats/simplified_radiosity_algorithm.py
View File

@ -67,7 +67,7 @@ class SimplifiedRadiosityAlgorithm:
i = (total_days + day - 1) * 24 + hour - 1
representative_building = self._city.buildings[0]
_global = representative_building.diffuse[cte.HOUR][i] / cte.WATTS_HOUR_TO_JULES
_beam = representative_building.beam[cte.HOUR][i] / cte.WATTS_HOUR_TO_JULES
_beam = representative_building.direct_normal[cte.HOUR][i] / cte.WATTS_HOUR_TO_JULES
content += f'{day} {month} {hour} {_global} {_beam}\n'
with open(file, 'w', encoding='utf-8') as file:
file.write(content)

4
hub/imports/energy_systems/montreal_future_energy_systems_parameters.py
View File

@ -82,8 +82,7 @@ class MontrealFutureEnergySystemParameters:
return _generic_energy_systems
@staticmethod
def _create_generation_systems(archetype_system):
def _create_generation_systems(self, archetype_system):
_generation_systems = []
archetype_generation_systems = archetype_system.generation_systems
if archetype_generation_systems is not None:
@ -107,6 +106,7 @@ class MontrealFutureEnergySystemParameters:
_generation_system.cell_temperature_coefficient = archetype_generation_system.cell_temperature_coefficient
_generation_system.width = archetype_generation_system.width
_generation_system.height = archetype_generation_system.height
_generation_system.tilt_angle = self._city.latitude
_generic_storage_system = None
if archetype_generation_system.energy_storage_systems is not None:
_generic_storage_system = ElectricalStorageSystem()

2
hub/imports/results/simplified_radiosity_algorithm.py
View File

@ -34,7 +34,7 @@ class SimplifiedRadiosityAlgorithm:
for key in self._results:
_irradiance = {}
header_name = key.split(':')
result = [x for x in self._results[key]]
result = [x * cte.WATTS_HOUR_TO_JULES for x in self._results[key]]
city_object_name = header_name[1]
building = self._city.city_object(city_object_name)
surface_id = header_name[2]

8
hub/imports/weather/epw_weather_parameters.py
View File

@ -114,10 +114,14 @@ class EpwWeatherParameters:
for x in self._weather_values['global_horizontal_radiation_wh_m2']]
building.diffuse[cte.HOUR] = [x * cte.WATTS_HOUR_TO_JULES
for x in self._weather_values['diffuse_horizontal_radiation_wh_m2']]
building.beam[cte.HOUR] = [x * cte.WATTS_HOUR_TO_JULES
for x in self._weather_values['direct_normal_radiation_wh_m2']]
building.direct_normal[cte.HOUR] = [x * cte.WATTS_HOUR_TO_JULES
for x in self._weather_values['direct_normal_radiation_wh_m2']]
building.beam[cte.HOUR] = [building.global_horizontal[cte.HOUR][i] -
building.diffuse[cte.HOUR][i]
for i in range(len(building.global_horizontal[cte.HOUR]))]
building.cold_water_temperature[cte.HOUR] = wh().cold_water_temperature(building.external_temperature[cte.HOUR])
# create the monthly and yearly values out of the hourly
for building in self._city.buildings:
building.external_temperature[cte.MONTH] = \

9
main.py
View File

@ -1,7 +1,6 @@
from scripts.geojson_creator import process_geojson
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
@ -18,8 +17,11 @@ from scripts.costs.cost import Cost
from scripts.costs.constants import SKIN_RETROFIT_AND_SYSTEM_RETROFIT_AND_PV, SYSTEM_RETROFIT_AND_PV
import hub.helpers.constants as cte
from hub.exports.exports_factory import ExportsFactory
import csv
from scripts.solar_angles import CitySolarAngles
from scripts.radiation_tilted import RadiationTilted
# Specify the GeoJSON file path
# geojson_file = process_geojson(x=-73.5953602192335, y=45.492414530022515, diff=0.001)
geojson_file = process_geojson(x=-73.5681295982132, y=45.49218262677643, diff=0.001)
file_path = (Path(__file__).parent / 'input_files' / 'output_buildings.geojson')
# Specify the output path for the PDF file
output_path = (Path(__file__).parent / 'out_files').resolve()
@ -30,9 +32,8 @@ city = GeometryFactory('geojson',
year_of_construction_field='year_of_construction',
function_field='function',
function_to_hub=Dictionaries().montreal_function_to_hub_function).city
# Enrich city data
# # Enrich city data
ConstructionFactory('nrcan', city).enrich()
UsageFactory('nrcan', city).enrich()
WeatherFactory('epw', city).enrich()
energy_plus_workflow(city)

21
scripts/pv_sizing_and_simulation.py Normal file
View File

@ -0,0 +1,21 @@
from scripts.radiation_tilted import RadiationTilted
from scripts.solar_angles import CitySolarAngles
import hub.helpers.constants as cte
class PVSizingSimulation:
def __init__(self, city, tilt_angle):
self.city = city
self.tilt_angle = tilt_angle
self.solar_angles = CitySolarAngles(file_name=self.city.name,
location_latitude=self.city.latitude,
location_longitude=self.city.longitude,
tilt_angle=self.tilt_angle,
surface_azimuth_angle=180,
standard_meridian=-75).calculate
self.enrich_radiation_data()
def enrich_radiation_data(self):
for building in self.city.buildings:
roof_horizontal = [x / cte.WATTS_HOUR_TO_JULES for x in building.roofs[0].global_irradiance[cte.HOUR]]
RadiationTilted(building, self.solar_angles, tilt_angle=self.tilt_angle, ghi=roof_horizontal).enrich()

109
scripts/radiation_tilted.py Normal file
View File

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

146
scripts/solar_angles.py Normal file
View File

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

2
tests/test_construction_factory.py
View File

@ -81,7 +81,7 @@ class TestConstructionFactory(TestCase):
self.assertEqual(len(building.external_temperature), 0, 'building external temperature is calculated')
self.assertEqual(len(building.global_horizontal), 0, 'building global horizontal is calculated')
self.assertEqual(len(building.diffuse), 0, 'building diffuse is calculated')
self.assertEqual(len(building.beam), 0, 'building beam is calculated')
self.assertEqual(len(building.direct_normal), 0, 'building beam is calculated')
self.assertIsNotNone(building.lower_corner, 'building lower corner is none')
self.assertEqual(len(building.sensors), 0, 'building sensors are assigned')
self.assertIsNotNone(building.internal_zones, 'no internal zones created')

2
tests/test_geometry_factory.py
View File

@ -52,7 +52,7 @@ class TestGeometryFactory(TestCase):
self.assertEqual(len(building.external_temperature), 0, 'building external temperature is calculated')
self.assertEqual(len(building.global_horizontal), 0, 'building global horizontal is calculated')
self.assertEqual(len(building.diffuse), 0, 'building diffuse is calculated')
self.assertEqual(len(building.beam), 0, 'building beam is calculated')
self.assertEqual(len(building.direct_normal), 0, 'building beam is calculated')
self.assertIsNotNone(building.lower_corner, 'building lower corner is none')
self.assertEqual(len(building.sensors), 0, 'building sensors are assigned')
self.assertIsNotNone(building.internal_zones, 'no internal zones created')

2
tests/test_usage_factory.py
View File

@ -44,7 +44,7 @@ class TestUsageFactory(TestCase):
self.assertEqual(len(building.external_temperature), 0, 'building external temperature is calculated')
self.assertEqual(len(building.global_horizontal), 0, 'building global horizontal is calculated')
self.assertEqual(len(building.diffuse), 0, 'building diffuse is calculated')
self.assertEqual(len(building.beam), 0, 'building beam is calculated')
self.assertEqual(len(building.direct_normal), 0, 'building beam is calculated')
self.assertIsNotNone(building.lower_corner, 'building lower corner is none')
self.assertEqual(len(building.sensors), 0, 'building sensors are assigned')
self.assertIsNotNone(building.internal_zones, 'no internal zones created')