wip using pvfactors engine simulation

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)

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
+    surface_azimuth = orientation["surface_azimuth"]
-    shaded_row_rotation = 90
+    axis_azimuth = orientation["axis_azimuth"]
-    shading_row_rotation = 90
-    surface_azimuth = 90  # east facing
-    axis_tilt = 0
-    axis_azimuth = 180
 
     gcr = np.divide(c["panel"]["dimensions"]["length"], pitch)
     gcr = min(1, gcr)
     logger.info(f"Ground coverage ratio: {gcr}")
 
-    # use pvfactors bifacial modelling package
+    pvrow_height = (
-    POA_data = pd.Series(dtype=object)
+        c["panel"]["dimensions"]["length"] * orientation["pvrow_height_ratio_to_length"]
-    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
 
-    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

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

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,
-    calculate_energy_production_vertical,
     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)