PVEngine needs to use timeseries albedo values (#98)

* Make sure engine uses all timeseries values of reflectivities

* Fix rho_mat in test

* Add test to catch issues with variable albedo

pvfactors/engine.py

@@ -195,46 +195,55 @@ def run_full_mode(self, fn_build_report=None):
         pvarray = self.pvarray
 
         # Get the irradiance modeling matrices
-        # shape = n_surfaces, n_timesteps
+        # shape = n_surfaces + 1, n_timesteps
         irradiance_mat, rho_mat, invrho_mat, _ = \
             self.irradiance.get_full_ts_modeling_vectors(pvarray)
 
         # --- Calculate view factors
-        # shape = n_surfaces, n_surfaces, n_timesteps
+        # shape = n_surfaces + 1, n_surfaces + 1, n_timesteps
         ts_vf_matrix = self.vf_calculator.build_ts_vf_matrix(pvarray)
         pvarray.ts_vf_matrix = ts_vf_matrix
         # Reshape for broadcasting and inverting
-        # shape = n_timesteps, n_surfaces, n_surfaces
+        # shape = n_timesteps, n_surfaces + 1, n_surfaces + 1
         ts_vf_matrix_reshaped = np.moveaxis(ts_vf_matrix, -1, 0)
 
         # --- Solve mathematical problem
-        # Build matrix of inverse reflectivities
-        # shape = n_surfaces, n_surfaces
-        invrho_mat = np.diag(invrho_mat[:, 0])
+        # shape of system = n_timesteps, n_surfaces + 1, n_surfaces + 1
+        shape_system = ts_vf_matrix_reshaped.shape
+        # Build 3d matrix of inverse reflectivities: diagonal for each ts
+        # shape = n_timesteps, n_surfaces + 1, n_surfaces + 1
+        invrho_ts_diag = np.zeros(shape_system)
+        for idx_surf in range(shape_system[1]):
+            invrho_ts_diag[:, idx_surf, idx_surf] = invrho_mat[idx_surf, :]
         # Subtract matrices: will rely on broadcasting
-        # shape = n_timesteps, n_surfaces, n_surfaces
-        a_mat = invrho_mat - ts_vf_matrix_reshaped
+        # shape = n_timesteps, n_surfaces + 1, n_surfaces + 1
+        a_mat = invrho_ts_diag - ts_vf_matrix_reshaped
         # Calculate inverse, requires specific shape
-        # shape = n_timesteps, n_surfaces, n_surfaces
+        # shape = n_timesteps, n_surfaces + 1, n_surfaces + 1
         inv_a_mat = np.linalg.inv(a_mat)
         # Use einstein sum to get final timeseries radiosities
-        # shape = n_surfaces, n_timesteps
+        # shape = n_surfaces + 1, n_timesteps
         q0 = np.einsum('ijk,ki->ji', inv_a_mat, irradiance_mat)
         # Calculate incident irradiance: will rely on broadcasting
-        # shape = n_surfaces, n_timesteps
-        qinc = np.dot(invrho_mat, q0)
+        # shape = n_surfaces + 1, n_timesteps
+        qinc = np.einsum('ijk,ki->ji', invrho_ts_diag, q0)
 
         # --- Derive other irradiance terms
         # shape = n_surfaces, n_timesteps
         isotropic_mat = ts_vf_matrix[:-1, -1, :] * irradiance_mat[-1, :]
         reflection_mat = qinc[:-1, :] - irradiance_mat[:-1, :] - isotropic_mat
 
         # --- Calculate AOI losses and absorbed irradiance
-        rho_mat = np.tile(rho_mat[:, 0], (rho_mat.shape[0], 1)).T
+        # Create tiled reflection matrix of
+        # shape = n_surfaces + 1, n_surfaces + 1, n_timestamps
+        rho_ts_tiled = np.moveaxis(np.tile(rho_mat.T, (shape_system[1], 1, 1)),
+                                   -1, 0)
+        # Get vf AOI matrix
         # shape [n_surfaces + 1, n_surfaces + 1, n_timestamps]
         vf_aoi_matrix = (self.vf_calculator
-                         .build_ts_vf_aoi_matrix(pvarray, rho_mat))
+                         .build_ts_vf_aoi_matrix(pvarray, rho_ts_tiled))
         pvarray.ts_vf_aoi_matrix = vf_aoi_matrix
+        # Get absorbed irradiance matrix
         # shape [n_surfaces, n_timestamps]
         irradiance_abs_mat = (
             self.irradiance.get_summed_components(pvarray, absorbed=True))

pvfactors/tests/test_engine.py

@@ -624,3 +624,42 @@ def test_engine_w_faoi_fn_in_irradiance_vfcalcs(params, pvmodule_canadian):
     np.testing.assert_almost_equal(
         pvarray.ts_pvrows[1].back.get_param_weighted('qabs'),
         [114.2143503])
+
+
+def test_engine_variable_albedo(params, df_inputs_clearsky_8760):
+    """Run PV engine calcs with variable albedo"""
+
+    irradiance_model = HybridPerezOrdered()
+    pvarray = OrderedPVArray.init_from_dict(params)
+    eng = PVEngine(pvarray, irradiance_model=irradiance_model)
+
+    # Manage timeseries inputs
+    df_inputs = df_inputs_clearsky_8760
+
+    # Get MET data
+    timestamps = df_inputs.index
+    dni = df_inputs.dni.values
+    dhi = df_inputs.dhi.values
+    solar_zenith = df_inputs.solar_zenith.values
+    solar_azimuth = df_inputs.solar_azimuth.values
+    surface_tilt = df_inputs.surface_tilt.values
+    surface_azimuth = df_inputs.surface_azimuth.values
+    albedo = np.linspace(0, 1, num=8760)
+
+    # Fit engine
+    eng.fit(timestamps, dni, dhi, solar_zenith, solar_azimuth,
+            surface_tilt, surface_azimuth, albedo)
+
+    # Run timestep
+    pvarray = eng.run_full_mode(fn_build_report=lambda pvarray: pvarray)
+    # Check the bifacial gain values
+    pvrow = pvarray.ts_pvrows[1]
+    bfg = (np.nansum(pvrow.back.get_param_weighted('qinc'))
+           / np.nansum(pvrow.front.get_param_weighted('qinc'))) * 100.
+    bfg_after_aoi = (np.nansum(pvrow.back.get_param_weighted('qabs'))
+                     / np.nansum(pvrow.front.get_param_weighted('qabs'))
+                     ) * 100.
+    expected_bfg = 16.441664
+    expected_bfg_after_aoi = 16.109509
+    np.testing.assert_allclose(bfg, expected_bfg)
+    np.testing.assert_allclose(bfg_after_aoi, expected_bfg_after_aoi)

pvfactors/tests/test_viewfactors/test_calculator.py

@@ -86,9 +86,9 @@ def test_vfcalculator_no_aoi_functions(params):
     assert np.sum(vfcalculator.vf_matrix) > 0
     # Compute reflectivity
     rho_mat = np.tile([0.03] * (pvarray.n_ts_surfaces + 1),
-                      (pvarray.n_ts_surfaces + 1, 1)).T
+                      (1, pvarray.n_ts_surfaces + 1, 1)).T
     assert rho_mat.shape == (pvarray.n_ts_surfaces + 1,
-                             pvarray.n_ts_surfaces + 1)
+                             pvarray.n_ts_surfaces + 1, 1)
     vf_aoi_matrix = vfcalculator.build_ts_vf_aoi_matrix(pvarray, rho_mat)
 
     # Check that correct size

pvfactors/viewfactors/calculator.py

@@ -136,7 +136,8 @@ def build_ts_vf_aoi_matrix(self, pvarray, rho_mat):
             PV array whose timeseries view factor AOI matrix to calculate
         rho_mat : np.ndarray
             2D matrix of reflectivity values for all the surfaces in the
-            PV array + sky. Shape = [n_ts_surfaces + 1, n_ts_surfaces + 1]
+            PV array + sky.
+            Shape = [n_ts_surfaces + 1, n_ts_surfaces + 1, n_timestamps]
 
         Returns
         -------
@@ -161,9 +162,7 @@ def build_ts_vf_aoi_matrix(self, pvarray, rho_mat):
         if self.vf_aoi_methods is None:
             # The reflection losses will be considered all diffuse.
             faoi_diffuse = 1. - rho_mat
-            # use broadcasting
-            vf_aoi_matrix = faoi_diffuse * np.moveaxis(vf_aoi_matrix, -1, 0)
-            vf_aoi_matrix = np.moveaxis(vf_aoi_matrix, 0, -1)
+            vf_aoi_matrix = faoi_diffuse * vf_aoi_matrix
         else:
             # Calculate vf_aoi between pvrow and ground surfaces
             self.vf_aoi_methods.vf_aoi_pvrow_to_gnd(ts_pvrows, ts_ground,