From 2504b4e5fc7d74d8eb10c2c2165beb2f3996256e Mon Sep 17 00:00:00 2001 From: Lucas Tan Date: Sun, 6 Apr 2025 16:08:44 +0800 Subject: [PATCH] wip using pvfactors engine simulation --- Utilities/Optimisation.py | 6 +- Utilities/Shading.py | 171 ++++++++++++++------------------------ config.yml | 12 ++- main.py | 7 +- 4 files changed, 81 insertions(+), 115 deletions(-) diff --git a/Utilities/Optimisation.py b/Utilities/Optimisation.py index b5fdce3..c313985 100644 --- a/Utilities/Optimisation.py +++ b/Utilities/Optimisation.py @@ -1,4 +1,4 @@ -from Utilities.Shading import calculate_energy_production_vertical +from Utilities.Shading import calculate_energy_production from scipy.optimize import minimize import logging @@ -19,7 +19,7 @@ def optimise_vertical_panel_pitch(c): pitch += c["panel"]["dimensions"]["thickness"] c["array"]["spacing"] = pitch logging.info(f"Optimizing with pitch: {pitch}m") - vertical_energy, _ = calculate_energy_production_vertical(c) + vertical_energy, _ = calculate_energy_production(c, "vertical") total_energy_yield = vertical_energy.sum() logger.info(f"Total energy yield for pitch {pitch}m: {total_energy_yield}kWh") return -total_energy_yield @@ -36,5 +36,5 @@ def optimise_vertical_panel_pitch(c): optimal_pitch = result.x[0] c["array"]["spacing"] = optimal_pitch logger.info(f"Optimal pitch found: {optimal_pitch}m") - vetical_energy, no_of_panels = calculate_energy_production_vertical(c) + vetical_energy, no_of_panels = calculate_energy_production(c, "vertical") return (optimal_pitch, vetical_energy, no_of_panels) diff --git a/Utilities/Shading.py b/Utilities/Shading.py index 5ac719b..d16f2f1 100644 --- a/Utilities/Shading.py +++ b/Utilities/Shading.py @@ -5,7 +5,9 @@ import math import matplotlib.pyplot as pl import pvlib -from pvlib.bifacial.pvfactors import pvfactors_timeseries +from pvfactors.geometry import OrderedPVArray +from pvfactors.engine import PVEngine +import matplotlib.animation as animation from Utilities.Processes import calculate_no_of_panels, calculate_required_system_size @@ -28,7 +30,7 @@ def define_grid_layout(c, panel_tilt): ) # calculate pitch - pitch = c["array"]["spacing"] + c["panel"]["dimensions"]["thickness"] + pitch = c["array"]["spacing"] # calculate minimum pitch if we don't want panel overlap at all min_pitch = c["panel"]["dimensions"]["length"] * math.cos( panel_tilt / 180 * math.pi @@ -136,9 +138,12 @@ def sanity_check_minimum_pitch(c): return solar_positions -def calculate_energy_production_vertical(c): +def calculate_energy_production(c, orientation): + orientation = c["array"]["orientation"][orientation] c = calculate_required_system_size(c) - panel_coordinates, no_of_panels = define_grid_layout(c, panel_tilt=90) + panel_coordinates, no_of_panels = define_grid_layout( + c, panel_tilt=orientation["panel_tilt"] + ) solar_positions, clearsky_data = get_solar_data(c) # the first row is always not shaded so exclude @@ -150,43 +155,66 @@ def calculate_energy_production_vertical(c): collector_width = c["panel"]["dimensions"]["length"] # calculate delta between unique y coordinates of panels to get pitch pitch = np.unique(panel_coordinates["y"])[1] - np.unique(panel_coordinates["y"])[0] - surface_to_axis_offset = 0 - shaded_row_rotation = 90 - shading_row_rotation = 90 - surface_azimuth = 90 # east facing - axis_tilt = 0 - axis_azimuth = 180 + surface_azimuth = orientation["surface_azimuth"] + axis_azimuth = orientation["axis_azimuth"] gcr = np.divide(c["panel"]["dimensions"]["length"], pitch) gcr = min(1, gcr) logger.info(f"Ground coverage ratio: {gcr}") - # use pvfactors bifacial modelling package - POA_data = pd.Series(dtype=object) - for row in range(0, 3): - result = pvfactors_timeseries( - solar_zenith=solar_positions["apparent_zenith"], - solar_azimuth=solar_positions["azimuth"], - surface_azimuth=surface_azimuth, - surface_tilt=90, - axis_azimuth=axis_azimuth, - timestamps=solar_positions.index, - dni=clearsky_data["dni"], - dhi=clearsky_data["dhi"], - gcr=gcr, - pvrow_height=c["panel"]["dimensions"]["length"], - pvrow_width=c["panel"]["dimensions"]["width"] * no_of_panels_in_row, - albedo=0.2, - n_pvrows=3, - index_observed_pvrow=row, - ) - # set negative values to 0 - poa_front = result[2].clip(lower=0) - poa_rear = result[3].clip(lower=0) - poa_global = poa_front + poa_rear * c["panel"]["bifaciality"] - POA_data.at[row] = poa_global + pvrow_height = ( + c["panel"]["dimensions"]["length"] * orientation["pvrow_height_ratio_to_length"] + ) - total_hourly_irradiance = POA_data.at[0] + POA_data.at[2] + (POA_data.at[1] * 20) + pvarray_parameters = { + "n_pvrows": 3, # number of pv rows + "pvrow_height": pvrow_height, # height of pvrows (measured at center / torque tube) + "pvrow_width": c["panel"]["dimensions"]["length"], # width of pvrows + "axis_azimuth": axis_azimuth, # azimuth angle of rotation axis + "gcr": gcr, # ground coverage ratio + } + + pvarray = OrderedPVArray.init_from_dict(pvarray_parameters) + engine = PVEngine(pvarray) + + inputs = pd.DataFrame( + { + "dni": clearsky_data["dni"], + "dhi": clearsky_data["dhi"], + "solar_zenith": solar_positions["zenith"], + "solar_azimuth": solar_positions["azimuth"], + "surface_tilt": np.repeat(orientation["panel_tilt"], len(solar_positions)), + "surface_azimuth": np.repeat(surface_azimuth, len(solar_positions)), + } + ) + inputs.index = clearsky_data.index + inputs.index.name = "index" + + engine.fit( + inputs.index, + inputs.dni, + inputs.dhi, + inputs.solar_zenith, + inputs.solar_azimuth, + inputs.surface_tilt, + inputs.surface_azimuth, + albedo=0.2, + ) + + pvarray = engine.run_full_mode(fn_build_report=lambda pvarray: pvarray) + + f, ax = pl.subplots(figsize=(10, 3)) + + def update(frame): + ax.clear() + pvarray.plot_at_idx(frame, ax, with_surface_index=True) + ax.set_title(inputs.index[frame]) + return ax + + ani = animation.FuncAnimation( + f, update, frames=len(inputs.index), interval=100, repeat=True + ) + pl.show() gamma_pdc = c["panel"]["temperature_coefficient"] temp_cell = c["panel"]["nominal_operating_cell_temperature"] p_row = no_of_panels_in_row * c["panel"]["peak_power"] @@ -217,76 +245,3 @@ def calculate_energy_production_vertical(c): logger.info(f"Total energy yield calculated: {total_energy} kWh") return total_hourly_energy, no_of_panels - - -def calculate_energy_production_horizontal(c): - c["array"]["system_size"] = ( - c["array"]["system_size"] * c["array"]["horizontal_max_capacity"] - ) - panel_coordinates, no_of_panels = define_grid_layout(c, panel_tilt=0) - solar_positions, clearsky_data = get_solar_data(c) - - # the first row is always not shaded so exclude - no_of_rows = np.unique(panel_coordinates["y"]).shape[0] - no_of_shaded_rows = no_of_rows - 1 - - collector_width = c["panel"]["dimensions"]["length"] - # calculate delta between unique y coordinates of panels to get pitch - pitch = np.unique(panel_coordinates["y"])[1] - np.unique(panel_coordinates["y"])[0] - surface_to_axis_offset = 0 - shaded_row_rotation = 0 - shading_row_rotation = 0 - axis_tilt = 0 - axis_azimuth = 270 # south facing surface - - projected_solar_zenith = pvlib.shading.projected_solar_zenith_angle( - solar_zenith=solar_positions["apparent_zenith"], - solar_azimuth=solar_positions["azimuth"], - axis_azimuth=axis_azimuth, - axis_tilt=axis_tilt, - ) - - shaded_fraction = pvlib.shading.shaded_fraction1d( - solar_zenith=projected_solar_zenith, - solar_azimuth=solar_positions["azimuth"], - axis_azimuth=axis_azimuth, - shaded_row_rotation=shaded_row_rotation, - shading_row_rotation=shading_row_rotation, - collector_width=collector_width, - pitch=pitch, - surface_to_axis_offset=surface_to_axis_offset, - axis_tilt=axis_tilt, - ) - shaded_fraction = shaded_fraction * no_of_shaded_rows / no_of_rows - - poa = pvlib.irradiance.get_total_irradiance( - surface_tilt=0, - surface_azimuth=180, - solar_zenith=projected_solar_zenith, - solar_azimuth=solar_positions["azimuth"], - dni=clearsky_data["dni"], - ghi=clearsky_data["ghi"], - dhi=clearsky_data["dhi"], - surface_type="urban", - ) - poa = poa.dropna(subset=["poa_global"]) - - effective_front = poa["poa_global"] * (1 - shaded_fraction) - - system_size = c["panel"]["peak_power"] * no_of_panels - - pdc0 = system_size - gamma_pdc = c["panel"]["temperature_coefficient"] - temp_cell = c["panel"]["nominal_operating_cell_temperature"] - pdc = pvlib.pvsystem.pvwatts_dc( - pdc0=pdc0, - gamma_pdc=gamma_pdc, - temp_cell=temp_cell, - g_poa_effective=effective_front, - ) - - total_hourly_energy = pdc * 15 / 60 / 1e3 # convert to kWh - total_energy = total_hourly_energy.sum() - logger.info(f"Total energy yield calculated: {total_energy} kWh") - - return total_hourly_energy, no_of_panels diff --git a/config.yml b/config.yml index 4585f63..a5ae1a5 100644 --- a/config.yml +++ b/config.yml @@ -5,8 +5,18 @@ array: edge_setback: 1.8 # distance from the edge of the roof to the array roof_slope: 0 slope: 0 # degrees from horizontal (+ve means shaded row is higher than the row in front) - horizontal_max_capacity: 0.75 # scale down due to peak power demand limit of NOVA performance_ratio: 0.9 # ratio of actual energy output to the theoretical maximum + orientation: + vertical: + surface_azimuth: 90 # degrees from North (clockwise) + axis_azimuth: 180 # degrees from North (clockwise) + panel_tilt: 90 + pvrow_height_ratio_to_length: 0.5 + horizontal: + surface_azimuth: 180 # degrees from North (clockwise) + axis_azimuth: 270 # degrees from North (clockwise) + panel_tilt: 0 + pvrow_height_ratio_to_length: 0 simulation_date_time: start: 2025-03-30 00:00 # start date and time in ISO 8601 format diff --git a/main.py b/main.py index 979b9da..6aca102 100644 --- a/main.py +++ b/main.py @@ -5,8 +5,7 @@ import numpy as np import matplotlib.pyplot as pl import matplotlib.dates as mdates from Utilities.Shading import ( - calculate_energy_production_horizontal, - calculate_energy_production_vertical, + calculate_energy_production, sanity_check_minimum_pitch, ) from Utilities.Optimisation import optimise_vertical_panel_pitch @@ -41,7 +40,9 @@ logger.debug(f"Vertical Energy Production: {vertical_energy.sum()}") logger.debug("Number of panels: %d", no_of_panels_vertical) logger.debug(f"System size: {no_of_panels_vertical * c['panel']['peak_power']/1e3} kWp") -horizontal_energy, no_of_panels_horizontal = calculate_energy_production_horizontal(c) +horizontal_energy, no_of_panels_horizontal = calculate_energy_production( + c, "horizontal" +) logger.info("Energy production for horizontal panels calculated successfully.") logger.debug(f"Horizontal Energy Production: {horizontal_energy.sum()}") logger.debug("Number of panels: %d", no_of_panels_horizontal)