Refactor to separate out VR computation within pvlib.scaling.wvm to a new function pvlib.scaling._compute_vr.

Author: jranalli

pvlib/scaling.py

@@ -62,26 +62,9 @@ def wvm(clearsky_index, positions, cloud_speed, dt=None):
 
     # Added by Joe Ranalli (@jranalli), Penn State Hazleton, 2019
 
-    pos = np.array(positions)
-    dist = pdist(pos, 'euclidean')
     wavelet, tmscales = _compute_wavelet(clearsky_index, dt)
 
-    # Find effective length of position vector, 'dist' is full pairwise
-    n_pairs = len(dist)
-
-    def fn(x):
-        return np.abs((x ** 2 - x) / 2 - n_pairs)
-    n_dist = np.round(scipy.optimize.fmin(fn, np.sqrt(n_pairs), disp=False))
-
-    # Compute VR
-    A = cloud_speed / 2  # Resultant fit for A from [2]
-    vr = np.zeros(tmscales.shape)
-    for i, tmscale in enumerate(tmscales):
-        rho = np.exp(-1 / A * dist / tmscale)  # Eq 5 from [1]
-
-        # 2*rho is because rho_ij = rho_ji. +n_dist accounts for sum(rho_ii=1)
-        denominator = 2 * np.sum(rho) + n_dist
-        vr[i] = n_dist ** 2 / denominator  # Eq 6 of [1]
+    vr = _compute_vr(positions, cloud_speed, tmscales)
 
     # Scale each wavelet by VR (Eq 7 in [1])
     wavelet_smooth = np.zeros_like(wavelet)
@@ -101,6 +84,68 @@ def fn(x):
     return smoothed, wavelet, tmscales
 
 
+def _compute_vr(positions, cloud_speed, tmscales):
+    """
+    Compute the variability reduction factors for each wavelet mode for the
+    Wavelet Variability Model [1-3].
+
+    Parameters
+    ----------
+    positions : numeric
+        Array of coordinate distances as (x,y) pairs representing the
+        easting, northing of the site positions in meters [m]. Distributed
+        plants could be simulated by gridded points throughout the plant
+        footprint.
+
+    cloud_speed : numeric
+        Speed of cloud movement in meters per second [m/s].
+
+    tmscales: numeric
+        The timescales associated with the wavelets in seconds [s].
+
+    Returns
+    -------
+    vr : numeric
+        an array of variability reduction factors for each tmscale.
+
+    References
+    ----------
+    .. [1] M. Lave, J. Kleissl and J.S. Stein. A Wavelet-Based Variability
+       Model (WVM) for Solar PV Power Plants. IEEE Transactions on Sustainable
+       Energy, vol. 4, no. 2, pp. 501-509, 2013.
+
+    .. [2] M. Lave and J. Kleissl. Cloud speed impact on solar variability
+       scaling - Application to the wavelet variability model. Solar Energy,
+       vol. 91, pp. 11-21, 2013.
+
+    .. [3] Wavelet Variability Model - Matlab Code:
+       https://github.com/sandialabs/wvm
+    """
+
+    # Added by Joe Ranalli (@jranalli), Penn State Hazleton, 2021
+
+    pos = np.array(positions)
+    dist = pdist(pos, 'euclidean')
+
+    # Find effective length of position vector, 'dist' is full pairwise
+    n_pairs = len(dist)
+
+    def fn(x):
+        return np.abs((x ** 2 - x) / 2 - n_pairs)
+
+    n_dist = np.round(scipy.optimize.fmin(fn, np.sqrt(n_pairs), disp=False))
+    # Compute VR
+    A = cloud_speed / 2  # Resultant fit for A from [2]
+    vr = np.zeros(tmscales.shape)
+    for i, tmscale in enumerate(tmscales):
+        rho = np.exp(-1 / A * dist / tmscale)  # Eq 5 from [1]
+
+        # 2*rho is because rho_ij = rho_ji. +n_dist accounts for sum(rho_ii=1)
+        denominator = 2 * np.sum(rho) + n_dist
+        vr[i] = n_dist ** 2 / denominator  # Eq 6 of [1]
+    return vr
+
+
 def latlon_to_xy(coordinates):
     """
     Convert latitude and longitude in degrees to a coordinate system measured

pvlib/tests/test_scaling.py

@@ -92,6 +92,13 @@ def expect_cs_smooth():
     return np.array([1., 1., 1.05774, 0.94226, 1.])
 
 
+@pytest.fixture
+def expect_vr():
+    # Expected VR for expecttmscale
+    return np.array([3., 3., 3., 3., 3., 3., 2.9997844, 2.9708118, 2.6806291,
+                     2.0726611, 1.5653324, 1.2812714, 1.1389995])
+
+
 def test_latlon_to_xy_zero():
     coord = [0, 0]
     pos_e = [0, 0]
@@ -172,6 +179,11 @@ def test_compute_wavelet_dwttheory(clear_sky_index, time,
     assert_almost_equal(np.sum(wavelet, 0), csi_series)
 
 
+def test_compute_vr(positions, expect_tmscale, expect_vr):
+    vr = scaling._compute_vr(positions, cloud_speed, np.array(expect_tmscale))
+    assert_almost_equal(vr, expect_vr)
+
+
 def test_wvm_series(clear_sky_index, time, positions, expect_cs_smooth):
     csi_series = pd.Series(clear_sky_index, index=time)
     cs_sm, _, _ = scaling.wvm(csi_series, positions, cloud_speed)